1 Presentación

Este documento es una versión actualizada -en varios sentidos- de unos materiales de análisis de microarrays que escribimos allá por 2010, usando LateX y Sweave, junto con la profesora M. Carme Ruíz de Villa para la asignatura “Análisis de Microarrays” del posgrado de Bioinformática de la Universitat Oberta de Catalunya .

En estos años han cambiado muchas cosas. El posgrado es ahora un máster interuniversitario entre la UOC y la Universidad de Barcelona (UB). Muchos estadísticos y/o bioinformáticos -entre los que me incluyo- han cambiado de LateX/Sweave a RStudio/RMarkdown y Github. Mi compañera de escritura, Mamen, se ha jubilado, pero un nuevo compañero, Ricardo, se ha apuntado a la aventura de revivir estos materiales. Y aunque sigue sin estar claro si los microarrays también se jubilaran pronto, el campo de las ómicas ha crecido de forma explosiva. El “Análisis de datos ómicos,” título de este nuevo documento ya no se centra en los microarrays. También trata de RNA-seq, Proteómica, Metagenómica, Metabolómica o Epigenómica por citar sólo algunos de los más populares.

Es obvio que un curso de “Análisis de datos ómicos” no puede abarcar todas estas disciplinas, aunque también lo es que comparten muchos elementos comunes.

El plan de desarrollo para esta nueva edición, que podrá seguirse en el repositorio de Github (), tendrá dos fases:

En una primera, actualizaremos la versión original a RMarkdown, eliminando aspectos desfasados y actualizando el código R a la versión actual de R/Bioconductor.
En una posterior miraré de añadir capítulos que introduzcan nuevas tecnologías y muestren como analizar los datos que estas generan. Si logro engañar a algunos colegas para que escriban un capítulo también los añadiré a los autores.

Previsiblemente este documento, con sus errores e imperfecciones se quedará como “open source” para uso de quien lo desee sin tener que pasar por el doloroso y anquilosante proceso de convertirlo en un libro, que, en este caso, empezaría a quedar obsoleto al minuto de estar acabado.

1.1 Instalación de paquetes de R y código de los capítulos

En este documento se utiliza cierto número de paquetes de R. En vez de indicar que se instalen la primera vez que aparecen hemos creado un archivo de instalación -que convertiremos más adelante en un paquete- llamado install_me_1st,R. Nuestra recomendación es que antes de empezar a usar los materiales se ejecute dicho archivo para asegurarse de que todos los paquetes necesarios se encuentran instalados, preferentemente en la versión más actualizada.

El archivo se encuentra en el directorio codigo_R.

Este directorio contiene además los archivos de código R extraídos de los archivos fuente “Rmarkdown” (con la instrucción knitr::purl()) para facilitar su ejecución.

2 Introduccion a los microarrays y otras tecnologías omicas

2.1 Antecedentes históricos

La biología molecular ha estado interesada desde sus comienzos en poder determinar el nivel de expresión de los genes integrados en el genoma humano. Para ello dispone desde hace años de múltiples técnicas para medir estos niveles, tales como el Northern blot (a nivel de ARN) o el Western blot (a nivel de proteína).

Especial interés ha tenido siempre sobre las otras biomoléculas el estudio de los niveles del ARN transcrito. Hace unos años se acuñó el término de transcriptómica para referirse al estudio del ARN en cada una de sus formas. La transcriptómica ha contribuido a tener una mejor comprensión de la patogénesis de las enfermedades.

En el año 2009 Wang y colaboradores (Wang, Gerstein, and Snyder (2009)) definieron el transcriptoma como el conjunto completo de transcritos en una célula y su cantidad en un momento determinado del desarrollo o en una determinada condición fisiológica. Como objetivos de la transcriptómica se han definido tres principales:

Catalogar todas las especies de transcritos, incluyendo ARNm, ARN no codificante y pequeños ARNs.
Determinar la estructura transcripcional de los genes, en términos de sus puntos 5’ de inicio y 3’ de finalización, las modificaciones pos-transcripcionales y los patrones de splicing.
Cuantificar los diferentes niveles de expresión de cada tránscrito durante el desarrollo y bajo diferentes condiciones.

El interés de la transcriptómica no solamente se ha centrado en el desarrollo de nuevas tecnologías que mejoran su estudio, sino también en el desarrollo de nuevos métodos de extraer la gran cantidad de información que se genera con estas nuevas técnicas.

2.2 Los microarrays

Una de las técnicas que revolucionó el estudio del transcriptoma fueron los microarrays.

Un microarray es un artefacto para la realización de experimentos que permite estudiar simultáneamente múltiples unidades, que que representan los genes, proteínas o metabolitos, sobre un sustrato sólido de cristal, plástico o sílice, y expuestos a la acción de las moléculas diana cuya expresión se desea analizar (ver Figura 2.1.

Figure 2.1: Imagen de un microarray

La utilización de los microarrays en la última década ha generado inmensas cantidades de datos útiles para el estudio y el desarrollo de enfermedades, dando a conocer múltiples mapas de expresión génica, encontrar biomarcadores o construir firmas génicas para determinadas enfermedades.

Lo que caracteriza a los nuevos métodos utilizados para estudiar el transcriptoma no es lo que pueden medir, sino la cantidad de mediciones simultáneas que pueden realizar. Mientras que hasta hace apenas una década se estudiaban los genes uno a uno en profundidad, a partir del uso de estas nuevas tecnologías se pueden estudiar muchísimos genes a la vez, pero en contrapartida con mucho menos detalle y más ruido.

Los microarrays, hoy una metodología bien consolidada, han sido cruciales para concebir una nueva manera de estudiar el transcriptoma, especialmente en el campo de la expresión génica. Después de ellos han venido otras técnicas que tienen en común con ellos el “alto rendimiento” es decir la capacidad para medir muchas variables –cientos o miles– a la vez. Entre estas técnicas podemos destacar los arrays de SNPs, los microRNAs, la metilación y especialmente la secuenciación de nueva generación.

2.2.1 Aplicaciones de los microarrays

La tecnología de los microarrays se ha aplicado a una inmensa variedad de problemas, desde el estudio de enfermedades como el cáncer o la esclerosis múltiple al de los ritmos circadianos de las frutas. Entre otros temas los microarrays se han aplicado a:

Estudio de genes que se expresan diferencialmente entre varias condiciones (sanos vs enfermos, mutantes vs salvajes, tratados vs no tratados) (Callow et al. (2000), Hedenfalk et al. (2002), Kupin (2005)).
Clasificación molecular en enfermedades complejas (Golub et al. (1999), Barrier et al. (2007)).
Identificación de genes característicos de una patología (firma o “signature”) (Dyrskjøt et al. (2002)).
Predicción de respuesta a un tratamiento (Lin, Devakumar, and Kibbe (2006)).
y una gran variedad de otros temas

2.3 Cómo funcionan los microarrays

En términos generales los microarrays funcionan mediante la hibridación de una sonda específica (“probe”) y una molécula diana (“target”). La hibridación que ha tenido lugar se detecta mediante fluorescencia y se visualiza con la ayuda de un escáner. Los niveles de fluorescencia detectados reflejan la cantidad de moléculas diana presentes en la muestra problema.

Existen diferentes tipos de microarrays según la naturaleza del “target” que se está estudiando. Podemos encontrar microarrays de proteínas, de tejidos, de ADN o de ARN (también llamados de expresión). En este curso nos centraremos en los microarrays de expresión génica aunque al final del capítulo se hace mención de algunos de los otros tipos de microarrays que existen.

Una clasificación habitual de los microarrays es por el número de muestras que se hibridan simultáneamente. Distinguimos:

Microarrays de dos colores o “spotted arrays”
Microarrays de un color o arrays de oligonucleótidos.

2.3.1 Microarrays de dos colores

Estos arrays aparecieron a mediados de los años 90 (Schena (1999)) y están basados en la hibridación competitiva de dos muestras, cada una de las cuales ha sido marcada con un marcador fluorescente diferente (normalmente Cy3 y Cy5, verde y rojo respectivamente).

En el array se imprimen las sondas (una sonda$\simeq$ un gen) cuyas secuencias se obtienen de información almacenada en bases de datos de secuencias como GenBank, dbEST, etc.

Figure 2.2: En los microrrays de dos colores las sondas son sintetizadas in vitro y depositadas directamente sobre una superficie de cristal

El ARN de las muestras problemas es extraído y posteriormente marcado con los marcadores fluorescentes, según a que grupo a comparar pertenezcan. Las muestras marcadas se mezclan y se hibridan sobre el array. La hibridación consiste en combinar la muestra (target) y el microarray (sondas) y dejarlos un tiempo en una cámara de hibridación a una temperatura y agitación determinadas. Las sondas que tengan secuencias complementarias en las muestras, se hibridarán con ellas y quedaran fuertemente adheridas. Pasadas unas horas se lava el microarray para eliminar los “targets” que no se hayan hibridado. Después de la hibridación el array se ilumina con un láser que provoca que el marcador fluorescente emita fluorescencia de uno u otro color generando dos imágenes que se superpondrán para su análisis conjunto. La cantidad de fluorescencia generada es proporcional a la cantidad de ARNm presente en la muestra problema. El resultado final es un valor que representa el nivel de expresión de una muestra respecto a la otra por lo que se le denomina expresión _ relativa_.

2.3.2 Microarrays de un color

Como su propio nombre indica, en estos arrays las muestras están marcadas únicamente con un marcador fluorescente. En cada array solamente se hibrida una muestra, por lo que no se da la hibridación competitiva como pasaba en los arrays de dos colores. El valor que se obtiene después de iluminar el array con el láser es una medida numérica que se obtiene directamente del escáner, es decir no está referida al valor de otra muestra por lo que recibe el nombre de expresión absoluta.

2.3.2.1 Fabricación de los arrays de oligonucleótidos

La empresa Affymetrix (Santa Clara, California) es la casa comercial líder en la manufacturación y venta de este tipo de microarrays. Affymetrix sintetiza las sondas directamente sobre el chip mediante un proceso llamado fotolitografía. Este proceso consiste en la adición cíclica de los cuatro nucleótidos (adenina, timina, citosina y guanina) sobre la superficie rígida, donde existen ancladas unas especies químicas reactivas que se protegen y desprotegen para añadir el nucleótido deseado, mediante ciclos de luz y oscuridad. Así se consigue la síntesis de oligonucleótidos de unos 25-mer (“mer”=bases) de longitud. En la figura (ref?)(fig:c01affymetryxsynthesis) se esquematiza el proceso.

Figure 2.3: Esquema del proceso de fabricación de arrays de oligonucleótidos

Cada oligonucleótido sintetizado de 25-mer (es decir 25 bases) se denomina sonda (“probe”). Alrededor de 40x107 de estas sondas se agrupan en una celda llamada “probe cell.” Las “probe cell” estan organizadas en parejas, “probe pairs,” donde se combinan un “Perfect Match” (PM) -cuya secuencia de 25-mer coincide perfectamente con una parte del gen-y un “Mismatch” (MM) -idéntico al PM, escepto en el nucléotido central que no coincide (ver la figura here)

De 11 a 20 “probe pairs” se agrupan para formar un “probe set.” Diferentes “probe sets” se distribuyen a lo largo del array, y son las que definen a qué gen se unen y las que definen los niveles de expresión detectados.

El Mismatch (MM) se creo originalmente para que proporcionara una medida de la unión inespecífica (es decir lo que se hibridara con esta secuencia “parecida” pero distinta a la original se consideraría ruído). Aunque algunos algoritmos originales de Affymetrix todavía lo utilizan, hoy en día ha caído en desuso, muchos algoritmos modernos como RMA –que se explicará más adelnate– ya no lo utilizan y en las versiones de arrays posteriores a 2005 ya no se incluyen.

Estructura de un grupo de sondas (`probe sets`) de un array antiguo de Affymetrix en el que se muestra los distintos términos definidos: `probe-set`, `probe-pair`, `Perfect-Match`, `Mismatch`

Figure 2.4: Estructura de un grupo de sondas (probe sets) de un array antiguo de Affymetrix en el que se muestra los distintos términos definidos: probe-set, probe-pair, Perfect-Match, Mismatch

2.3.3 Realización de un experimento de microarrays

Tanto en los arrays de un color como en los de dos colores el material de partida es el ARN. Es muy importante controlar la calidad y la cantidad del ARN, ya que puede influir directamente en la calidad final de los resultados. Para ello se utiliza un equipo llamado “Bioanalyzer” (Agilent Technologies, Santa Clara, California) que proporciona entre otros parámetros un valor llamado “RNA Integrity Number (``RIN”)’’ que sirve para decidir si la calidad del ARN es suficientemente buena como para que valga la pena usarlo para hibridar un microarray. Este valor varia entre 0 y 10 y en general suele exigirse un valor mínimo de por ejemplo 7 u 8.

Una vez decicido que la calidad del ARN de la muestra es aceptable la cantidad inicial de ARN es amplificada unas 25 veces utilizando enzimas y cebadores específicos.

Posteriormente el ARN es marcado con una molécula fluorescente (ficoeritrina en este caso) y fragmentado en trozos más pequeños que se puedan unir a las “probes” fijadas en el array, mediante el proceso de hibridación. Al igual que en los arrays de dos colores, el array hibridado es iluminado en un escáner mediante un láser, para que la ficoeritrina emita luz fluorescente. La fluorescencia registrada de cada “probe set,” es directamente proporcional a la cantidad de ARN presente en la muestra inicial de cada gen. La figura here es una representación esquemática del uso de microarrays de un color.

Figure 2.5: Esquema del funcionamiento de un microarrays de Affymetrix o de un sólo color

2.3.4 Como se mide la expresión

Los microarrays permiten cuantificar la expresión de los genes a través de la intensidad de la fluorescencia que es capturada por los escáneres. Las imágenes se convierten en valores por un proceso que no es objeto de discusión en este curso, pero que puede ser considerado relativamente fiable y estable (véase, por ejemplo Schena (1999)).

Cada tecnología genera diferentes tipos de imágenes y éstas generan diferentes valores que deben ser tratados adecuadamente para proporcionar alguna clase de estimaciones de una misma variable: la expresión génica.

Tal como se ha indicado anteriormente una de las principales diferencias entre arrays de uno y dos colores es que estos últimos se basan en la hibridación competitiva de las dos muestras mientras que los de un color sólo miden cuanta muestra se hibrida con las sondas del chip. En consecuencia en los arrays de dos colores se mide cuánto se expresa un gen en una muestra respecto a la otra lo que tiene un sentido biológico en términos de “sobre–expresión” (por ejemplo, si un gen se expresa dos veces más en una condición que en la otra puede deducirse que está activado o sobre–expresado) o “sub–regulación” (si el gen se expresa la mitad en una condición que en otra puede deucirse que está inhibido o regulado negativamente). En los arrays de un color en cambio tan s’olo se mide cuánto se expresa un gen en una escala que carece de sentido biológico.

2.3.4.1 Medición de la expresión relativa en arrays de dos colores

Cuando la imagen obtenida en un microarray de dos colores es analizada (“cuantitativamente”), en cada spot se generan algunos valores (ver figura here). Aunque esto depende del software utilizado, básicamente consiste en (i) Medidas de señales, Rojo ($R$) o Verde ($G$) para cada canal , (ii) Medidas de background, $R_b$ , $G_b$, que intentan proporcionar una medida de la fluorescencia no debida a hibridación, y (iii) algunas medidas de calidad para el spot.

Estas cantidades pueden utilizarse para proporcionar medidas sencillas del ratio de expresión: \[ M=\frac{R}{G}, \] o el background corregido del ratio de expresión: \[ M=\frac{R-R_b}{G-G_b}. \]

Es habitual la utilización del logaritmo en base 2 de esta cantidad como el resultado final de expresión relativa. Esto se debe principalemte a dos motivos: por un lado, los datos de expresión se aproximan mejor con una distribución log–normal, y por otro lado, la utilización de logaritmos simetriza las diferencias haciendo más fácil la interpretación.

Figure 2.6: Esquema del tratamiento de las imágenes para su cuantificación

La tabla here muestra de forma simplificada el aspecto que puede tener una matriz de expresión relativa en en la que para estudiar 6 muestras apareadas, tres de tejidos sanos y tres tumorales, se ha hibridado cada tumor contra su control “normal” dando lugar a una matriz de expresión de 5 genes y 3 columnas de expresiones relativas.

kableExtra::kable(tabexpresRel)

genes	logTum1Norm1	logTum2Norm2	logTum3Norm3
Gene 1	0.46	0.80	1.51
Gene 1	-0.90	0.06	-3.20
Gene 1	0.15	0.04	0.09
Gene 1	0.60	1.06	1.35
Gene 1	-0.45	-1.03	-0.79

2.3.4.2 Medición de la expresión absoluta en arrays de un color

Los arrays de de Affymetrix representan cada gen como un conjunto de sondas cada una de las cuales se corresponde con una fragmento corto de un gen. De hecho, tal como se ha dicho, no se trata de sondas sinó de parejas de sondas formadas por un “Perfect Match” que corresponde a la cadena de ADN original y un “Mismatch” en las que ha cambiado su nucleotido central.

La idea subyacente en esta aproximación es que cualquiera que hibrida con la sonda de mismatch no debería representar una expresión real sino que representa lo que se denomina “background” o señal de fondo.

En los inicios de ésta tecnología la compañía Affymetrix sugirió combinar ambas medidas en lo que puede verse como una medida de expresión corregida para la señal de fondo. La fórmula utilizada ha evolucionado, pero, una estimación sencilla de las primeras versiones es: \[ Avg.diff=\frac{1}{|A|}\sum_{j \in A}(PM_j-MM_j), \] donde $A$ es el conjunto de pares de sondas cuyas intensidades no se desvian más de tres veces de la desviación estandar de la principal intensidad entre todas las sondas.

La tabla here muestra de forma simplificada el aspecto que puede tener una matriz de expresión relativa en en la que para estudiar 6 muestras apareadas, tres de tejidos sanos y tres tumorales, se ha hibridado cada tumor y cada control “normal” en un array separado dando lugar a una matriz de expresión de 5 filas (una por gen) y 6 columnas de expresiones absolutas.

kableExtra::kable(tablexpresAbs)

genes	logTum1	logTum2	logTum3	logNorm1	logNorm2	logNorm3
Gene 1	5.10	6.90	6.60	6.40	7.80	4.30
Gene 2	6.97	8.74	7.89	8.03	9.70	5.63
Gene 3	4.44	6.89	6.41	6.02	7.47	4.08
Gene 4	6.43	8.13	8.56	7.14	7.63	4.81
Gene 5	8.18	10.44	13.29	7.85	9.60	5.29

La mayor diferencia entre estas dos maneras de medir la expresión no está en la fórmula específica, que ha evolucionado en ambos casos, sino en el hecho que, mientras que en los chips de Affymetrix se tiene un único valor de expresión para cada condición, en los arrays de dos colores se trabaja con una medida de la expresión relativa entre dos condiciones. Aunque Affymetrix permite estimaciones más precisas, les expresiones relativas tienen una mejor interpretación intuitiva.

2.4 Bioinformática de Microarrays

El crecimiento en el uso de la experiencia en microarrays en la última década ha sido en paralelo por los desarrollos necesarios en la metodología –nuevos métodos para modelar y analizar los datos se requieren a menudo– y la bioinformática –nuevas herramientas necesarias para implementar los métodos, así como el almacenamiento, el acceso o la organización de la mayor parte creciente de datos disponibles. Esto nos lleva a considerar dos aspectos muy importantes relacionados con el análisis de microarrays de datos:

¿Qué software está disponible para analizar los datos de microarrays?
¿Qué sistemas de base de datos están disponibles para almacenar y manipular los datos de microarrays tanto a nivel local como global?

Este tema puede ser considerado complementario, pero necesario para poner en práctica los puntos discutidos en el documento de manera que una breve presentación de los sistemas de software y base de datos se presentan a continuación.

2.4.1 Software para el análisis de datos de microarrays

Supongamos que un estadista quiere involucrarse en el análisis de datos de microarrays y después de leer un poco este entiende lo que hay que hacer. Una pregunta obvia es “¿Qué herramienta debo usar?”Como la mayoría de los profesionales en el campo que está familiarizado con varios paquetes y probablemente tiene algunas preferencias.

Después de algunas búsquedas en Google, es obvio que haya varias posibilidades

Para utilizar los paquetes estadísticos estándar –SPSS o SAS– y analizar los datos que deben haber sido procesados o exportados del texto.
Para usar una de las muchas herramientas disponibles gratuitamente, ya sea web o de base local.
Para contar con extensiones específicas de microarrays para el análisis de datos tales como el Proyecto Bioconductor.
Para comprar uno de los programas comerciales existentes.

Como es habitual cada opción tiene aspectos positivos y negativos. El uso de paquetes estadísticos estándar –SPSS o SAS– tiene la curva de aprendizaje más corta, pero no permite hacer la mayoría de los pre–pasos de procesamiento, tales como la normalización o resumen, por lo que debe combinarse con otro software. Además, si uno desea hacer un ANOVA o K–means están bien, pero si lo que uno quiere hacer es aplicar métodos específicos, tales como SAM o ajustes locales–FDR serán insuficientes. Algunos paquetes de estadística como S+ o ha desarrollado extensiones de gran alcance para el análisis de datos de microarrays.

2.4.1.1 Software libre o de código abierto

Existen numerosas herramientas gratuitas para el análisis de Microarrays. Algunas de ellas, además, llevan la posibilidad (libertad) de modificar el programa por que incluyen el codigo fuente de programa (por tanto, se las llamas de “software de código abierto” o de “software libre,” ver: http://www.gnu.org/philosophy/free-sw.html ). Y otras, en cambio, ponen restricciones a como permiten usar el programa, restringen la libertad de modificar el código fuente del programa, y por tanto, aunque sean gratuitas, se las llama de software propietario, no-libre o cerrado, por contraposicion al software anterior.

Cualquiera de estas herramientas gratuitas puede conllevar un esfuerzo considerable en aprender a utilizarlas, pues tienen diferentes grados de madurez, y de facilidad para que nuevos usuarios las utilicen. Asimismo, puede pasar que esta herramientas no sigan estándares comunes de organización de los datos o flujos de trabajo, lo que el aprendizaje de uno no suele ayudar al aprendizaje de otro. Y cabe la posibilidad de que, como herramientas gratuitas o demasiado jóvenes, presenten una mayor tasa de errores que lo deseado. En cualquier caso, a menudo pueden resultar útiles para un análisis exploratorio de nuestros datos de microarrays, o para el mundo de la enseñanza, en especial si hay una comunidad de usuarios lo suficientemente grande de esa herramienta para garantizar que no hay demasiados problemas básicos con su uso. Pero si se desea usarlos repetidamente para la realización de estudios de media a alta complejidad, pueden resultar ser insuficiente, ya sea porque carecen de métodos, porque son ineficientes o simplemente porque no tienen la capacidad de programación para automatizar tareas repetitivas.

Algunos listados de esta herramientas son:

A pesar de estas críticas, los programas libres especialmente (más allá de los únicamente “gratuitos”), pueden ser una forma suave para introducirse a los análisis de datos de microarrays. Para guiar a un usuario inexperto comentamos brevemente algunas de nuestras herramientas libres favoritas.

_BRB array tools_ es un complemento de Excel-ins que combina R, C y Java para hacer los cálculos y utiliza Excel para interactuar con el usuario – lo que significa que sólo está disponible para usuarios de Windows. Es proporcionada por The Biometrics Research Branch del Instituto Nacional del Cáncer (EE.UU.). éste se complementa con tutoriales y una base de datos y de estudios para estudios reales preparados para ser utilizados con estos. Suele ser muy atractivo a primera vista, especialmente cuando se utiliza con sus propios ejemplos. Sin embargo la creación de un nuevo análisis desde el principio no es una tarea fácil y lo peor es que tiende a chocar en una dura forma con mensajes de Visual Básic crípticos, especialmente si se utilizan en los ordenadores con las versiones no inglesas de Windows.
_TM4_ es un conjunto de cuatro programas de libre escritos en Java y se ejecutados en sistemas Linux y Windows desarrollados por el instituto TIGR (ahora J. Craig Venter). Aunque un poco viejo y relativamente sesgado hacia los arrays de dos colores, para el cual fue desarrollado originalmente, es muy robusto (se bloquea mucho menos que BRB) y ofrece capacidades de análisis, no sólo (MeV), sino también el análisis de imágenes (Spotfinder), la normalización separado (MIDAS) y un sistema de base de datos (MADAM) para almacenar los experimentos.
Un grave inconveniente de las herramientas anteriores es su sesgo histórico hacia los microarrays de dos colores lo que implica que se pierden (a partir de comienzos de 2008) los métodos de preprocesamiento importante como RMA. Una buena –fácil de usar– alternativa para los primeros pasos de control de calidad y pre–procesamiento de los chips de Affymetrix es ofrecido por la misma compañía. Se llama Consola de expresión y se puede descargar desde la web de Affymetrix después de la inscripción gratuita.
_The http://babelomics.bioinfo.cipf.es/_ es un conjunto integrado de herramientas para el análisis de datos de microarrays disponible en la web. Babelomics ha sido diseñado para proporcionar una intuitiva interfaz basada en Web que ofrece diversas opciones de análisis de la etapa inicial de pre–procesamiento (normalización de Affymetrix y experimentos de microarrays de dos colores y otras opciones de pre–procesamiento), hasta el paso final de la elaboración de perfiles funcionales de la experiencia (con gene Ontology, pathways, PubMed resúmenes, etc), que incluyen diferentes posibilidades de agrupación, de predicción de genes de selección de clase y serie-comparativo de gestión de la hibridación genómica.

{r c01GEPAS, fig.align='center', fig.cap="Un mapa de las funcionalidades de "Babelomics" organizados como en una línea de metro. Un usuario normalmente debe comenzar en alguna parte de la izquierda del mapa y finalizar en algún lugar de la derecha", echo=FALSE} knitr::include_graphics("epsimages/c01GEPAS.png")

2.4.1.2 El proyecto bioconductor

Una de las opciones para el análisis de los datos mencionados anteriormente es la combinación de un software estándar, tales como Matlab, Mathematica o R, con librerías específicas diseñadas para el análisis de microarrays. Aunque existen algunas extensiones de Matlab para el análisis de microarrays es con R, que esta complementariedad ha alcanzado dimensiones inesperadas. El proyecto Bioconductor (www.Bioconductor.org) comenzó en 2001 como un proyecto de código abierto y libre desarrollo de software para el análisis y comprensión de datos genómicos. Su gran éxito ha hecho crecer a partir de poco más de una docena de paquetes a cientos de ellos. Casi todas las técnicas disponibles en el análisis de microarrays tienen su propio paquete, y con frecuencia hay varios de ellos.

La gran potencia de este proyecto implica también algunos de sus inconvenientes: en primer lugar, al ser un proyecto de código abierto significa que los desarrolladores contribuyen con sus programas a “tal cual.” Aunque hay sistemas de control para evitar la no– ejecución del código, es más difícil de garantizar (a excepción de la honestidad de los desarrolladores) que se ejecuta como se indica. El poder del Bioconductor también se basa en la flexibilidad de la lengua R. Es muy difícil para los usuarios que no son competentes R hacer un uso eficiente de estas bibliotecas.

A pesar de estas aparentes dificultades, Bioconductor es la herramienta elegida por muchos estadísticos y la razón principal es que, cuando uno ha sido capaz de sentirse cómodo con ella, su poder es difícil de igualar. Las instalaciones de la programación de R, hacen posible automatizar el análisis, así como la generación de informes, por lo que es la opción de elección cuando se realizan tareas repetitivas.

2.4.1.3 El software propietario

Hay muchas herramientas comerciales disponibles para el análisis de microarrays de datos. Estos van desde pequeños programas específicos de un tipo de datos en paquetes de software grandes, como Partek Genomics Suite que es una solución completa optimizada para los cálculos eficientes y rápidos, así como para la mayoría de los datos genómicos existentes. Software comercial de microarrays tiene los pros y los contras del tradicional software comercial: Puede ser bueno, pero es caro y puede que no sea suficientemente flexible para un usuario experto que desea introducir sus propios métodos en el análisis.

2.4.2 Bases de datos de Microarrays

La diversidad de formatos y tipos de microarrays de experimentos ha hecho difícil que un formato de base de datos cualquiera se haya impuesto y no hay sistema de base de datos que se haya convertido en el “dorado–estándar.”

En efecto, existe como estado algún tipo de acuerdo sobre la información mínima sobre un experimento de microarrays que debe ser almacenada (the MIAME standard (www.mged.org/Workgroups/MIAME/miame.html) es un acrónimo de esto), pero como si fuera un tema político que el acuerdo haya sido tan corto que es más simbólico que útil.

Se pueden distinguir dos niveles en los que los sistemas de bases de datos se han desarrollado.

Sistemas de bases de datos locales El análisis de los datos de microarrays pasa por una serie de pasos en los diferentes tipos de datos, imágenes, archivos binarios, archivos de texto tienen que ser procesados. Se requiere que tener almacenados en un lugar fácilmente accesible. Algunos sistemas como BASE(base.thep.lu.se/) o caArray (caarray.nci.nih.gov/) son soluciones de gran alcance para el almacenamiento de datos y experimentos, pero su uso está lejos de ser tan extendido como el de las herramientas de software de análisis.
Repositorios públicos de arrays La comunidad biológica se ha comprometido, desde el inicio de los microarrays, que los datos de los experimentos publicados deben hacerse públicos. Esto ha creado la necesidad de repositorios de microarrays pública donde cualquier usuario puede almacenar sus datos en una forma adecuada. Al mismo tiempo se ha hecho una impresionante cantidad de datos disponibles para un nuevo análisis para cualquiera que desee hacerlo, que ofrece una riqueza sin precedentes de oportunidades cuyo poder está empezando a mostrar. Una lista de las colecciones de datos pública disponible en http://www.nslij-genetics.org/microarray/data.html.

2.5 Extensiones (1): Otros tipos de microarrays

Este curso se centrará, en torno al tipo más popular de microarrays: los microarrays de expresión de ARN, diseñados para estudiar la expresión génica basada en la información sobre la cantidad de ADN que se transcribe el ARNm.

La disponibilidad de las tecnologías del genoma ha permitido desarrollar otros tipos de microarrays. Por “otros,” se puede decir que se basan en microarrays de ADN para estudiar otros problemas de microarrays de expresión o que dependen de otras sustancias como las proteínas o los hidratos de carbono. Una descripción completa de cada tipo, su uso, objetivos y análisis de los datos está fuera del alcance de este trabajo. Sin embargo, para dar un ejemplo de las similitudes y diferencias entre los microarrays de expresión y las tecnologías relacionadas hacemos una breve reseña de los problemas que requieren de estas tecnologías alternativas y dar una breve descripción de uno de ellos, los arrays de genotipasdo o “de SNPs.”

2.5.1 Otros microarrays de ADN

Uno de los ejes principales de la genómica funcional es la comprensión y la curación de la enfermedad. Se sabe que muchas alteraciones genéticas subyacen anomalías y/o enfermedades. Por ejemplo:

Mutaciones puntuales –cambio de una o más bases– puede dar lugar a proteínas alteradas o cambios en el nivel de expresión.
La pérdida de copias de genes puede reducir el nivel de expresión. Estos cambios están relacionados con la supresión de tumores.
Ganancia de copias del gen puede aumentar el nivel de expresión y están con la activación de oncogenes relacionados.
Metilación o de–metilación de los promotores de genes, respectivamente, pueden disminuir o aumentar el nivel de expresión. Estos también están relacionados con los oncogenes supresores de tumores.

Los diferentes tipos de microarrays se adaptan a estudiar las manifestaciones y los efectos de estas alteraciones. Los puntos planteados anteriormente pueden ser estudiados con (i) o SNP y (ii) hibridación comparativa del genoma o microarrays de ADN CGH y otros como o arrays de promotores.

2.5.1.1 Arrays de genotipado (SNPs)

El polimorfismo de un nucleótido es una forma de mutación puntual que consiste en las variaciones de pares de bases individuales que se encuentran dispersos al azar por todo el genoma. Miles de polimorfismos de un nucleótido han sido –y siguen siendo– identificados como parte de los proyectos de secuenciación del genoma. SNPs han sido altamente conservados durante la evolución y dentro de una población. Debido a esta medida de conservación el mapa de SNPs sirve como un excelente marcador genotípico para la investigación.

Los arrays SNP son un tipo de chips de ADN para detectar polimorfismos dentro de las poblaciones. Estos trabajan bajo los mismos principios básicos que arrays de expresión, pero cada sonda está diseñada para detectar las diferentes variaciones de polimorfismos de nucleótido único para cada SNP conocidos.

Figure 2.7: Utilización de arrays de genotipado

Los arrays SNP tienen muchas aplicaciones. Entre ellos se puede destacar:

Estudios basados en el linaje familiar El ADN de los familiares afectados con una enfermedad en particular puede ser comparado con el ADN de los miembros de la misma familia que no tienen esta condición. Estos estudios permiten identificar las diferencias genéticas que pueden estar asociados con la enfermedad.
Estudios de asociación a nivel poblacional consiste en determinar las diferencias en las frecuencias de SNP en los individuos afectados y no afectados en una población. El objetivo es identificar SNPs en particular o combinaciones de SNP que difieren entre los dos grupos y por lo tanto, asociados con la enfermedad. Estos estudios requieren un gran número de muestras para representar adecuadamente a la población. Este es uno de los mejores - aplicación conocida de las matrices de SNPs que ilustra cómo se puede ayudar en la identificación de genes relacionados con enfermedades complejas.
Cambios del número de copias SNPs pueden ser utilizados como etiquetas para las regiones con variabilidad del número de copias –una variante del número de copias (CNV)– es un segmento de ADN que es de 1 KB o más grande y está presente en un número variable de copias en comparación con un genoma de referencia’’. La identificación de los cambios de número de copias es útil para detectar tanto las aberraciones cromosómicas como la pérdida del número de copias neutrales de heterocigosidad (LOH), eventos que son característicos de muchos tipos de cáncer.

2.5.2 Otros tipos de microarrays: Proteínas o carbohidratos

Existe un amplio consenso sobre el hecho de que la información obtenida de los microarrays de ADN no es suficiente para alcanzar una completa comprensión de los procesos celulares la mayoría de los cuales son controlados por las proteínas que interactúan con otras moléculas como los hidratos de carbono a menudo implicados en importantes mecanismos biológicos como la interacción de patógenos con el huésped, el desarrollo o la inflamación. Las microarrays de proteínas (i) y carbohidratos (ii) son dos ejemplos de la extensión del uso de estas herramientas para el análisis de alto rendimiento de diferentes tipos de moléculas. Microarrays de tejidos (iii) son un tipo diferente de extensión en la que la resta no es distintas variantes de un solo tipo de molécula, sino de un tipo de tejidos.

2.6 Extensiones (2): NGS: Next(Now) Generation Sequencing

A finales de la primera década del siglo XXI parece ser que la revolución de los microarrays está llegando a su fin ((ledford:2008?), frantz:2005) y una nueva tecnología está inundando las revistas y congresos científicos. Se trata de la ultrasecuenciación también llamada “Next Generation Sequencing” o “Now Generation Sequencing.”

Básicamente la ultrasecuenciación permite hacer lo mismo que la secuenciación tradicional o “Sanger,” es decir obtener la secuencia de una cadena de ADN o ARN que se ha preparado previamente para ello mediante creación de un cierto número de copias o “librerías.” La diferencia, una vez más, está en la cantidad de secuencias que es posible obtener. Mientras que la secuenciación tradicional suele poder producir XXX pares de bases por día la capacidad de los nuevos métodos es órdenes de magnitud superior lo que significa el acceso rápido y económico a la capacidad de secuenciar genomas de eucariotas en un tiempo breve -semanas o días en 2012- a un coste ínfimo de secuenciar.

Este gran incremento en la capacidad de producir secuencias ha representado, de nuevo, un cambio de paradigma en la resolución de muchos problemas científicos. Hoy es posible considerar una aproximación directa al estudio de múltiples problemas (“Si quieres saber como va secuéncialo”) y la secuenciación se aplica a un gran número de campos desde la metagenómica que estudia a nivel genómico la diversidad de los ecosistemas bacterianos como lagos o intestinos, al análisis de variantes estructurales en exomas o genomas -y su posible asociación con enfermedades herediatarias como el autismo o el cáncer de mama. Una variante de la secuenciación de ADN ha sido el RNA-seq que, secuenciando ADN complementario permite abordar un estudio mucho más fino de la expresión génica de la que se pueda hacer con microarrays, dado que permite cuantificar directamente el producto de la expresión -podemos secuenciar el ARN y contar los transcritos para saber cuanta expresión se ha dado- a la vez que es posible obtener información acerca de secuencias no identificadas previamente -por lo que no se habían podido poner en los microarrays.

En esta sección se presenta una visión general de la ultrasecuenciación, las tecnologías que funcionan en 2011, los problemas bioinformáticos y computacionales que aparecen y como se resuelven y algunos de los problemas que pueden abordarse con estas técnicas con un sesgo hacia los aspectos bioinformáticos y de análisis de datos que no dejan de ser el objeto de este curso.

2.6.1 Visión global de las tecnologías de secuenciación

El ADN no puede ser secuenciado de golpe, ha de ser fragmentado suficientemente, secuenciado a trozos y finalmente reensamblado. La característica básica presente, tanto en el paradigma tradicional (Sanger) como en ultrasecuenciación es que, previa fragmentación de ADN a secuenciar, cada fragmento deberá ser amplificado o clonado. Debido a este paso previo de fragmentación y amplificación se puede acuñar a las tecnologías de secuenciación en general con el término shotgun sequencing o shotgun cloning. En cierto modo mayor amplificación implica mayor precisión de lectura. Dado un fragmento concreto a secuenciar y un cliclo de proceso, se intervendrá químicamente en su colonia o grupo de cadenas clones asociada (polony) de una manera particular dependiendo de la tecnología concreta. La primera diferencia entre la tecnología convencional (Sanger) y las NGS es qe la primera necesita que el resultado de la intervención sea diferente para cada cadena de la polony, mientras que las segundas todo lo contrario. Al final de este capítulo entenderemos que la precisión del resultado depende de la eficiencia conseguida en uno u otro sentido.

Desde principios de los 90, la producción de secuenciación de ADN ha sido llevada a cabo casi exclusivamente mediante implantaciones automáticas, y con escasa capacidad de procesamiento paralelo, de la bioquímica de Sanger. [4, 5].

2.6.2 De la secuenciación Sanger a la ultrasecuenciación

El método Sanger evolucionó desde su versión más artesanal hasta su versión automatizada actual. La clave principal del método de Sanger artesanal fue el uso de didesoxinucletidos trifosfato (ddNTPs) como terminadores de cadena, motivo por el cual también se hace referencia a éste con el término . Para la lectura de un solo fragmento, deben realizarse por separado, cuatro mezclas de reacción. Cada mezcla de reacción contiene los cuatro nucleótidos trifosfato (dATP, dCTP, dTTP . dGTP), ADN polimerasa I, un cebador o primer marcado radiactivamente (permite el acoplamiento de la polimerasa e inicio de la reacción en ese punto) y un nucleótido didesoxi ddNTP (A,C,G, o T), a una concentración baja. El nucleótido didesoxi utilizado competirá con su homólogo por incorporarse a los clones del fragmento dado, produciendo la terminación de la síntesis en el momento y lugar donde se incorpora, ya que al carecer de grupo 3’-OH no es posible el enlace fosfodiester con el siguiente nucleótido. Por este sistema, en cada mezcla de reacción se producen una serie de moléculas de ADN de nueva síntesis de diferente longitud que terminan todas en el mismo tipo de nucleótido y marcadas todas radiactivamente por el extremo 5’ (todas contienen en el extremo 5’ el primer utilizado).

Las secuecnias de ADN de nueva síntesis obtenidos en cada mezcla de reacción se separan por tamaños mediante electrofresis en geles verticales de acrilamida muy finos (0.5 mm de espesor) y de gran longitud (cerca de 50 cm) que permiten distinguir secuencias de ADN que se diferencian en un solo nucleótido. Los productos de cada una de las cuatro mezclas de reacción se insertan en cuatro carriles diferentes del gel.

Una vez terminada la electrofresis, el gel se pone en contacto con una película fotográfica de autorradiografía. La aparición de una banda en una posición concreta de la autorradiografía en uno de los cuatro carriles nos indica que en ese punto de la secuencia del ADN de nueva síntesis (complementario al ADN molde) está la base correspondiente al nucleótido didesoxi utilizado en la mezcla de reacción correspondiente.

Teniendo en cuenta que el ADN de nueva síntesis crece en la dirección 5’ -3’, si comenzamos a leer el gen por las secuencias de menos tamaño (extremo 5’) y avanzamos aumentando el tamaño de las secuencias (hacia 3’), obtendremos la secuencia completa del ADN de nueva síntesis en la dirección 5’ -3’ del fragmento dado.

La principal diferencia entre el método enzimático de terminación de cadena y el método automático de secuenciación radica, en primer lugar, en el tipo de marcaje. En el método automático en vez de radiactividad se utiliza fluorescencia y lo habitual es realizar cuatro mezclas de reacción, cada una con nucleótido trifosfato (dTTP) marcado con un fluorocromo distinto. Este sistema permite automatizar el proceso de manera que es posible leer al mismo tiempo los ADNs de nueva síntesis producto de las cuatro mezclas de reacción.

Después de tres décadas de perfeccionamiento gradual, la bioquímica Sanger puede ser aplicada actualmente para alcanzar longitudes de lectura de hasta 1000 pares de bases (“bp”) y precisiones crudas por base de hasta un 99.999% aunque con un coste de 0.5 $ por kilobase.

2.6.3 Tecnologías de ultrasecuenciación

Las distintas estrategias NGS pueden ser agrupadas dentro de varias categorías. Nosotros centramos nuestro interés en el desarrollo posterior únicamente en la categoría . El resto de categorías corresponden a tecnologías en desarrollo y quedarán muy brevemente introducidas en la siguiente lista:

Microchip-based electrophoretic sequencing: se trata de un intento de incremento de automatización y paralelismo de los métodos tradicionales; todavía en desarrollo, idealmente, se pretende integrar todos los aspectos del procesado de muestras, lo que reduciría sustancialmente el tiempo de proceso y el consumo de reactivos sin renunciar a la precisión del método tradicional.
Real-time observation of single molecules: una de las aproximaciones en este sentido viene de la mano de Pacific Biosciences; debido a la restricción de las reacciones de iluminación acopladas a la actividad de la polimerasa a un volumen del orden del zeptolitro ($10^{-21}$L), están consiguiendo que la incorporación de los nucleótidos pueda ser observada con muy bajo ruido; incluye un potencial de longitud de lectura mayor que el de Sanger con una capacidad enorme de paralelismo.
Sequencing by hybridization : el concepto básico de esta aproximación es que las diferencias de hibridación de ciertos fragmentos de ADN etiquetados contra un array de sondas de oligonucleótidos pueden ser usadas para identificar con precisión posiciones de variabilidad (variantes).
Cyclic-array sequencing:
- secuenciación 454 (usada en los 454 Genome Sequencers, Roche Applied Science; Basel).
- tecnología Solexa (usada en el Illuminia (San Diego) Genome Analyzer).
- la plataforma SOLiD (Applied Biosystems; Foster City, CA, USA).
- Polonator (Dover/Harvard).
- tecnología HeliScope Single Molecule Sequencer (Helicos; Cambridge, MA, USA). }

Aunque las plataformas incluidas en nuestra categoría de interés (cyclic-array sequencing) difieren en su bioquímica, el flujo de trabajo es conceptualmente similar (Fig. (ref?)(fig:sangevsngs) a). Hay varias aproximaciones para la generación de los clusters de clones (polonies), pero todas ellas suponen que al final los derivados PCR de una única molécula de la librería acaben clusterizados en un soporte espacial. Entre ellas están las in situ polonies (sujeción de localización única), bridge PCR (sujeción a sustrato plano) y emulsion PCR (sujeción a gotas (beads) de tamaño en torno al micrómetro ($10^{-6}$m) o micron). El proceso de secuenciación (véase fig: (ref?)(fig:sangevsngs) b) consiste en general en alternar ciclos de irrigación enzimática con adquisición de datos basada en procesado de imagen.

Figure 2.8: Comparación entre la secuenciación Sanger y la ultrasecuenciación

3 El proceso de análisis de microarrays

Este capítulo es una corta transición entre la primera parte del curso, en la que se han presentado los conceptos y herramientas básicos y la segunda en donde se presentan por separado y con mayor detalle los métodos de análisis de datos de microarrays.

Su objetivo por tanto es ofrecer una visión de conjunto que sirva de guía (“roadmap”) para los capítulos siguientes de forma que sin perder el detalle de cada uno de ellos tengamos conciencia de en que punto del proceso general nos encontramos.

El capítulo se estructura en dos partes. En la primera se presentan brevemente algunos de los problemas que típicamente se puede querer estudiar con microarrays u otras técnicas similares de análisis de datos de alto rendimiento. A continuación se presenta lo que se ha llamdo aquí el proceso de análisis de microarrays. Finalmente se intoducen algunos casos reales que, a modo de ejemplo se utilizarán en los capítulos siguientes.

3.1 Tipos de estudios

Los microarrays y otras tecnologías de alto rendimiento se han aplicado a multitud de investigaciones, de tipus muy diversos que van desde estudio del cancer (Alizadeh et al. (2000), Golub et al. (1999), van ’t Veer et al. (2002)) al de la germinación y la maduración del tomate (Moore et al. (2002)). A pesar de ello no resulta complicado clasificar los estudios realizados en algunos de los grandes bloques que se describen a continuación. La clasificación está basada en el excelente texto de Simon y colegas (Simon et al. (2003)) y aunque se origina en problemas de microarrays se puede aplicar fácilmente a estudios de genómica o ultrasecuenciación.

3.1.1 Comparación de grupos o “Class comparison”

El objetivo de los estudios comparativos es determinar si los perfiles de expresión génica difieren entre grupos previamente identificados. También se conoce estos estudios como de selección de genes diferencialmente expresados y son, sin duda los más habituales. Los grupos pueden representar una gran variedad de condiciones, desde distintos tejidos a distintos tratamientos o múltiples combinaciones de factores experimentales.

El análisis de este tipo de experimentos, que se describe en el capítulo sobre selección de genes diferencialmente expresados utiliza herramientas estadísticas como las pruebas de comparación de grupos paramétricas (t de Student) o no (test de Mann-Whitney) o diversos métodos de análisis de la varianza.

Entre los ejemplos de la sección @ref{c4examples} los casos @ref{celltypes}, @ref{estrogen} o @ref{CCL4} hacen referencia a estudios comparativos.

3.1.2 Predicción de clase o “Class prediction”

La predicción de clase puede confundirse con la selección de genes en tanto que disponemos de clases predefinidas pero su objetivo es distinto, ya que no pretende simplemente buscar genes cuya expresión sea distinta entre dichos grupos sino genes que puedan ser utilizados para identificar a que clase pertenece un “nuevo” individuo dado cuya clase es “a priori” desconocida. El proceso de predicción suele empezar con una selección de genes informativos, que pueden ser, o no, los mismos que se obtendrían si aplicáramos los métodos del apartado anterior, seguida de la construcción de un modelo de predicción y, lo que es más importante, de la verificación o validación de dicho modelo con unos datos nuevos independientes de los utilizados para el desarrollo del modelo.

Aunque el interés de la predicción de clase es muy alto se trata de un procedimiento mucho más complejo y con más posibilidades de error que la simple selección de genes diferencialmente expresados.

Entre los ejemplos de la sección @ref{c4examples} el caso @ref{golub} trata de un problema de predicción, a la vez que uno de descubrimiento de clases.

3.1.3 Descubrimiento de clases o “Class discovery”

Un problema distinto a los descritos se presenta cuando no se conoce las clases en que se agrupan los individuos. En este caso de lo que se trata es de encontrar grupos entre los datos que permitan reunir a los individuos más parecidos entre si y distintos de los de los demás grupos. Los métodos estadísticos que se emplearan en estos casos se conocen como análisis de clusters y no son tan complejos como los de predicción de clase aunque algunos aspectos como por ejemplo la definición del número de grupos no resulta tampoco sencillo.

Entre los ejemplos de la sección tanto el caso golub, en parte, como el @ref{breastTum} tratan problemas de descubrimiento de clases.

Una curiosidad del campo de la estadística es que el término clasificación aparece usado de forma indistinta para referirse a problemas de predicicón de clase o de descubrimiento de clase.

3.1.4 Otros tipos de estudios

Una vez identificados los principales tipos de estudios quedan muchos que no coinciden plenamente con ninguno de ellos. Sin entrar en detalles podemos señalar los estudios de evolución a lo largo del tiempo (“time course”), los de significación biológica (“Gene Enrichment Analysis,” “Gene Set Enrichment Analysis,” …) , los que buscan relaciones entre los genes (“network analysis” o “pathway analysis”). De momento con conocer e identificar los tres grandes bloques mencionados resultará más que suficiente.

3.2 Algunos ejemplos concretos

Una de las dificultades con que se encuentra la persona que comienza en el análisis de datos de microarrays es de donde obtener ejemplos concretos con los que prácticar las técnicas que está aprendiendo.

No es difícil encontrar datos de microarrays en internet por lo que se han seleccionado algunos conjuntos de datos interesantes para utilizarlos de ejemplo a lo largo del curso. Algunos de éstos son “populares” en el sentido de que han sido utilizados en diversas ocasiones y por lo tanto se encuentran bien documentados. Otros se han escogido simplemente porque ilustran bien algunas de las ideas que se desea exponer o por su accesibilidad.

Todos los datos corresponden a investigaciones publicadas por lo que no se describen exhaustivamente sinó que se expone brevemente el origen y objetivos del trabajo –incluyendo su clasificación según los grupos definidos en la sección anterior– y las características perinentes para el análisis como el tipo de microarrays, los grupos –si los hay– o como acceder a los datos.

3.2.1 Estudio de procesos regulados por citoquinas

Este conjunto de datos, que se denominará “celltypes,” corresponde a un estudio realizado por Chelvarajan y sus colegas (Chelvarajan et al. (2006)) que analizaron el efecto de la estimulación con lipopolisacáridos en la regulación por parte de citoquinas de ciertos procesos biológicos relacionados con la inflamación.

Este estudio es del tipo “class comparison” es decir su principal objetivo es la obtención de genes diferencialmente expresados entre dos o más condiciones.

Los datos se encuentran disponibles en la base de datos pública GEne Expression Omnibus mantenida por el National Institute of Health (NIH), pero pueden descargarse de la página de materiales del curso para garantizar su disponibilidad.

3.2.2 Clasificación molecular de la leucemia

A finales de los años 90, Todd Golub y sus colaboradores (Golub et al. (1999)) realizaron uno de los estudios más populares hasta el momento con datos de microarrays. En él utilizaron microarrays de oligonucleótidos para 6817 genes humanos para mirar de encontrar una forma de distinguir (clasificar) tumores de pacientes con leucemia linfoblástica aguda (ALL) de aquellos que sufrían de leucemia mieloide aguda (AML). Además se interesaba por la posibilidad de descubrir subgrupos de forma que pudieran definirse variantes de cada una de estas patologías a nivel molecular.

La diversidad de objetivos del estudio lleva a clasificarlo tanto entre los del tipo de predicción de clase como entre los que buscan descubrir nuievas clases o grupos en los datos.

Los datos de este estudio se encuentran disponibles en la web del instituto Broad, en donde se llevó a cabo (http://www.broadinstitute.org). También se encuentra disponible un paquete de R denominado ALL que permite utilizarlos directamente usando R y Bioconductor.

3.2.3 Efecto del estrogeno y el tiempo de administración

Scholtens y colegas (Scholtens et al. (2004)) describen un estudio sobre el efecto de un tratamiento con estrógenos y del tiempo transcurrido desde el tratamiento en la expresión génica en pacientes de cáncer de mama. Los investigadores supusieron que los genes asociados con una respuesta temprana podrían considerarse dianas directas del tratamiento, mientras que los que tardaron más en hacerlo podrían considerarse objetivos secundarios correspondientes a dianas más alejadas en las vías metabólicas.

Estos datos han sido utilizados multitud de veces en los cursos de análisis de microarrays realizados por el proyecto Bioconductor y se encuentran disponibles en forma de paquete de R, el paquete estrogen. Una cracaterística importante de este paquetes es el hecho de que en vez de los datos procesados proporciona los datos “crudos”en forma de archivos .CEL de Affymetrix. Esto permite una mayor flexibilidad a la hora de reutilizarlos lo que explica su popularidad.

3.2.4 Efecto del CCL4 en la expresión génica

Holger y colegas de la empresa LGC Ltd. en Teddington, Inglaterra realizaron unos expeimentos con microarrays de dos colores en los que trataron hepatocitos de ratón con tetraclorido de carbono (CCL4) o con dimetilsulfóxido (DMSO). El tetraclorido de carbono fue ampliamente utilizado en productos de limpieza o refrigeración para el hogar hasta que se detectó que podía tener efectos tóxicos e incluso cancerígenos. El DMSO es un solvente similar, sin efectos tóxicos conocidos, que se utilizó como control negativo.

Los datos de este estudio no han sido publicados pero se encuentran disponibles en el paquete CCL4 de bioconductor. Su interés reside por un lado en que se trata de datos de microarrays de dos colores de la marca Agilent –en un momento en que la mayoría de estudios se realizan con datos de un color. Aparte de esto cabe resaltar el hecho de que el paquete incluye, de forma similar al anterior, los datos “crudos” en forma de archivos de tipo “Genepix” uno de los programas populares para escanear imágenes generadas con microarrays de dos colores.

Este estudio es también un estudio de comparación de clases cuyo objetivo principal es la selección de genes cuya expresión se asocia al tratamiento con CCL4 o DMSO.

3.2.5 Estudio de la mutación “swirl” en el desarrollo del pez cebra

3.2.6 Análisis de patrones en el ciclo celular

Los datos de este ejemplo denominado yeast son datos ya normalizados referidos a la expresión de los genes en distintos momentos del ciclo delular de la levadura e decir desde que concluye la división celular hasta que se inicia la siguiente.

Los datos puede descargarse desde la página del proyecto “Yeast Cell Cycle Project” (Proyecto de estudio del ciclo celular de la levadura) en la dirección:

http://genome-www.stanford.edu/cellcycle/data/rawdata/.

3.2.7 Recapitulación

La tabla resume la lista de ejemplos que se utilizan en este manual indicando el nombre con que nos referiremos de aquí en adelante a cada conjunto de datos as' como algunas de sus carcaterísticas.

La tabla siguiente resume los “datasets” utilizados en este manual. Aparte del nombre (arbitrario y “mnemotécnico”) se indica el tipo de microarrays, el número de muestras, y el tipo de problema para el que se utilizaron originalmente.

Nombre	Tipo (Modelo)	N. Muestras	Tipo de estudio
`celltypes`	Un color (Affy, Mouse 4302)	12	Comparativo
`golub`	Un color (Affy, HGU95A)	38	Clasificación
`estrogen`	Un color (Affy, HGU95A)	8	Comparativo
`CCL4`	Dos colores (Agilent, WG Rat Microarray)	16	Comparativo
`swirl`	Dos colores (“Custom”)	4	Comparativo
`breastTum`	Un color (Affy, HGU95A)	49	Clasificación

3.3 El proceso de análisis de microarrays

Una vez descritos los tipos de análisis y algunos ejemplos podemos pasar a describir el proceso de análisis de microarrays que se resume brevemente en la figura @ref{c04analysisProcess}.

El análisis de microarrays, como la mayoría de análisis debe proceder de forma ordenada y siguiendo el método científico:

La pregunta y su contexto nos servirán de guía para definir el Diseño experimental adecuado.
El experimento se deberá realizar siguiendo las pautas decididas en el Diseño experimental y los datos obtenidos, que solemos denominar datos “crudos” o “raw data,” deberán someterse a los Controles de calidad adecuados antes de continuar con su análisis.
Una vez decidido si la calidad de los datos es aceptable pasaremos a prepararlos para el análisis lo que puede incluir diversas formas de preprocesado, o transformaciones que a menudo se incluyen de forma general bajo el paraguas del término normalización, aunque, como veremos se trata de conceptos distintos.
Los datos normalizados se utilizarán para los análisis estadísticos que hayamos decidido realizar durante el diseño experimental.
Finalmente los resultados de los análisis serán la base para una interpretación biológica de los resultados del experimento.

Tal como ilustra la figura @ref(fig:c04analysisProcess>) el análisis de microarrays puede ser fácilmente visualizado como un proceso que empieza por una pregunta biológica y concluye con una interpretación de los resultados de los análisis que, de alguna forma, confiamos nos acerque un poco a la respuesta de la pregunta inicial.

Figure 3.1: Proceso de análisis de microarrays

El proceso descrito es básicamente una forma razonable de proceder en general. Los microarrays y otros datos genómicos son diferentes en su naturaleza de los datos clásicos alrededor de los que se han desarrollado la mayor parte de técnicas estadísticas. En consecuencia, en muchos casos ha sido necesario adaptar las técnicas existentes o desarrollar otras nuevas para adecuarse a las nuevas situaciones encontradas. Esto ha determinado que existan muchos métodos para cada una de los pasos descritos anteriormente lo que da lugar a una grandísima cantidad de posibilidades.

En la práctica lo que suele hacerse es optar por utilizar algunos de los métodos en los que hay un cierto consenso acerca de su calidad y utilidad para cada problema. Allison (Allison et al. (2006)) repasa los puntos principales de este consenso dando una lista de puntos a tener en cuenta en cualquier estudio que utilice microarrays. Imbeaud y Auffray (Imbeaud and Auffray (2005)) citan una lista de hasta 39 puntos que uno debe seguir en un experimento con microarrays para usar “buenas prácticas.”

Finalmente Zhu y otros (Zhu, Miecznikowski, and Halfon (2010)) utilizan un conjunto de arrays con valores de expresión conocidos para proponer los que, a su parecer, resultan los métodos más apropiados para cada etapa desde la corrección del backround hasta la selección de genes diferencialmente expresados. La figura 3.2 ilustra algunas de las opciones sugeridas por dichos autores.

Figure 3.2: Diseño de arrays

Los capítulos que siguen al presente proceden aproximadamente en el orden del proceso descrito en 3.1.

Se empieza por tratar los principios del diseño de experimentos.
A continuación se describen algunos métodos para el control de calidad, el preprocesado y la normalización de los datos.
Se sigue con los métodos de selección de genes, que tratan aproximaciones básicas como el t-test y otras más sofisticadas como los modelos lineles para microarrays.
La parte de selección de genes finaliza con una introducción a los métodos de análisis de la significación biológica de las listas de genes obtenidas de los procesos anteriores.
Un último bloque aborda los métods de clasificación supervisada y no supervisada, que han tenido y siguen teniendo gran aplicabilidad en el análisis de todo tipo de datos ómicos.

4 Diseño de experimentos de microarrays

En este capítulo se examinan unas componentes clave del análisis de microarrays el diseño de experimentos que, no tan solo es crucial para una buena recogida de información, sino que marca todo el proceso desde el preprocesado al análisis final.

4.1 Fuentes de variabilidad

Los datos genómicos son muy variables. La figura 4.1, tomada de Geschwind and Gregg (2002) ilustra algunas posibles fuentes de variabilidad, desde que se inicia el experimento hasta que se lee la información con el escáner. Sin tener que entrar en cómo influye cada una de estas fuents o cuan grande es su influencia lo que muestra esta figura es que este tipo de experimentos se enfrenta a múltiples fuentes de error, por lo que puede beneficiarse de un correcto diseño de experimentos.

Figure 4.1: Proceso de análisis de microarrays

4.2 Tipos de variabilidad

Habitualmente, en la mayor parte de situaciones experimentales, podemos distinguir entre variabilidad sistemática y variabilidad aleatoria.

La variabilidad sistemática es principalmente debida a procesos técnicos mientras que la variabilidad aleatoria es atribuible tanto a razones técnicas como biológicas. Podemos encontrar ejemplos de variabilidad sistemática en la extracción del ARN, marcaje o foto detección. La variabilidad aleatoria puede estar relacionada con muchos factores tales como la calidad del ADN o las características biológicas de las muestras.

La manera natural de manejar la variabilidad aleatoria es, por supuesto, la utilización de un diseño experimental apropiado y el apoyo para obtener conclusiones de unas herramientas estadísticas adecuadas. Los problemas relacionados con el c de los experimentos los discutiremos en esta sección y los relacionados con la aplicación de métodos estadísticos serán tratados en la sección métodos estadísticos.

Tradicionalmente, se estiman las correcciones de la variabilidad sistemática a partir de los datos, en lo que se llama genéricamente “calibrado.” En este contexto, hablaremos de “normalización,” que se tratará más adelante, en el capítulo dedicado al preprocesado de los datos.

4.3 Conceptos básicos de diseño de Experimentos

Empezaremos por definir conceptos que aparecen de forma usual cuando se plantea un diseño:

Unidad experimental: Entidad física a la que se aplica un tratamiento, de forma independiente al resto de unidades. En cada unidad experimental se pueden realizar una medida o varias medidas, en este caso distinguiremos entre unidades experimentales y unidades observacionales.
Factor: son las variables independientes que pueden influir en la variabilidad de la variable de interés.
- Factor tratamiento: es un factor del que interesa conocer su influencia en la respuesta.
- Factor bloque: es un factor en el que no se está interesado en conocer su influencia en la respuesta, pero se supone que esta existe y se quiere controlar para disminuir la variabilidad residual.
Niveles : cada uno de los resultados de un factor. según sean elegidos por el experimentador o elegidos al azar de una amplia población se denominan factores de efectos fijos o factores de efectos aleatorios.
Tratamiento : es una combinación especifica de los niveles de los factores en estudio. Son, por tanto, las condiciones experimentales que se desean comparar en el experimento. En un diseño con un único factor son los distintos niveles del factor y en un diseño con varios factores son las distintas combinaciones de niveles de los factores.
Tamaño del Experimento: es el número total de observaciones recogidas en el diseño.

4.4 Principios de diseño de experimentos

Al planificar un experimento hay tres principios básicos que se deben tener siempre en cuenta: La replicación, el control local o “bloqueo” y la aleatorización.

Aplicados correctamente dichos principios garantizan diseños experimentales eficientes.

4.4.1 Replicación

La replicación o repetición de un experimento de forma idéntica en un número determinado de unidades, es la que permite la realización de un posterior análisis estadístico.

La utilización de replicas es importante para incrementar la precisión, obtener suficiente potencia en los tests y como base para los procedimientos de inferencia. Normalmente, distinguimos dos tipos de replicas en el análisis de microarrays:

La replicación técnica se utiliza cuando estamos tratando réplicas del mismo material biológico. Puede ser tanto la replicación de spots en el mismo chip como la de diferentes alícuotas de la misma muestra hibridadas en diferentes microarrays.
La replicación biológica se da cuando se toman medidas sobre muestras independientes, normalmente sobre individuos diferentes.

La replicación técnica permite la estimación del error a nivel de medida, mientras que la replicación biológica permite estimar la variabilidad a nivel de población.

Figure 4.2: Replicas biológicas vs réplicas técnicas

4.4.1.1 Potencia y tamaño de la muestra

Sorprendentemente, los primeros experimentos de microarrays utilizaban pocas o ninguna replica biológica. La principal explicación para este hecho, además de la falta de conocimiento estadístico, fue el alto coste de cada microarray.

En pocos años, la necesidad de las réplicas ha llegado a ser indiscutible, y al mismo tiempo, el coste de los chips ha disminuido considerablemente. Actualmente, es normal la utilización de, al menos, tres a cinco replicas por condición experimental, aunque a este consenso se ha llegado más por razones empíricas que por la disponibilidad de modelos apropiados para el análisis de la potencia y tamaño de la muestra.

Recientemente, se ha producido una importante afluencia de artículos describiendo métodos para el análisis de la potencia y tamaño de la muestra. A pesar de su variedad, ningún método aparece como candidato claro para ser utilizado en situaciones prácticas. Esto se debe, probablemente, a la complejidad de los datos de microarrays, básicamente por que los genes no son independientes. Por tanto, estas estructuras de correlación existen en los datos, pero la mayor parte de las dependencias son desconocidas por lo que la estimación de estas estructuras es muy complicada.

Como indicaba Allison Allison et al. (2006), aunque no hay consenso sobre que procedimiento es mejor para determinar el tamaño de las muestras, sí que lo hay sobre la conveniencia de realizar el análisis de la potencia, y, por supuesto, sobre el hecho de que un mayor número de replicas generalmente proporcionan una mayor potencia.

4.4.1.2 Pooling

En el contexto de los microarrays, llamamos “pooling” a la combinación del RNA de diferentes casos en una única muestra.

Hay dos razones a favor de ello, una, es que, a veces, no hay suficiente ARN disponible y esta es la única forma de conseguir suficiente material para construir los arrays, otra, más controvertida, es la creencia que la variabilidad entre arrays puede reducirse por “pooling.” La justificación es que combinar muestras equivale a promediar expresiones, y como ya se conoce, el promedio es menos variable que los valores individuales. A pesar de la debilidad de este argumento, es verdad que en ciertas situaciones el pooling puede ser apropiado y, recientemente, muchos estadísticos han dedicado sus esfuerzos para tratar de responder la pregunta “to pool or not to pool” (Kerr and Churchill 2001). Por ejemplo, si la variabilidad biológica está altamente relacionada con el error en las medidas, y las muestras biológicas tienen un coste mínimo en comparación con el de los arrays, una apropiada estrategia de “pooling” puede ser claramente eficiente en costes.

De todos modos, el “pooling” no se debería usar en cualquier tipo de estudios. Si el objetivo es comparar expresiones medias (ver más adelante ” comparación de clases”), puede funcionar adecuadamente, pero se debería claramente evitar cuando el objetivo del experimento es construir predictores que se basen en características individuales.

4.4.2 Aleatorización

Se entiende por aleatorizar la asignación de todos los factores al azar a las unidades experimentales. Co ello se consigue disminuir el efecto de los factores no controlados por el experimentador en el diseño experimental y que podrían influir en los resultados.

Las ventajas de aleatorizar los factores no controlados son:

Transforma la variabilidad sistemática no planificada en variabilidad no planificada o ruido aleatorio. Dicho de otra forma, aleatorizar previene contra la introducción de sesgos en el experimento.
Evita la dependencia entre observaciones al aleatorizar los instantes de recogida muestral.
Valida muchos de los procedimientos estadísticos más comunes.

4.4.3 Bloqueo o control local

Hace referencia a dividir o particionar las unidades experimentales en grupos llamados bloques de modo que las observaciones realizadas en cada bloque se realicen bajo condiciones experimentales lo más parecidas posibles. A diferencia de lo que ocurre con los factores tratamiento, el experimentador no está interesado en investigar las posibles diferencias de la respuesta entre los niveles de los factores bloque.

Bloquear es una buena estrategia siempre y cuando sea posible dividir las unidades experimentales en grupos de unidades similares. La ventaja de bloquear un factor que se supone que tiene una clara influencia en la respuesta, pero en el que no se está interesado, es que convierte la variabilidad sistemática no planificada en variabilidad sistemática planificada.

4.5 Diseños experimentales para microarrays de dos colores

En arrays de dos colores, se aplican dos condiciones experimentales a cada array. Esto permite la estimación del efecto del array, como el efecto de bloqueo. En Affymetrix u otro array de un canal, cada condición se aplica a un chip separado, imposibilitando la estimación del efecto de los microarrays, lo cual, por otra parte, se considera que tiene una relación muy pequeña en el tratamiento de los efectos debido al proceso industrial utilizado para fabricar estos chips.

Como consecuencia de lo anterior, los experimentos que usan arrays de un canal se consideran “estándar,” por lo que se les pueden aplicar los conceptos y técnicas tradicionales de diseño de experimentos .

Los arrays de dos canales presentan una situación más complicada. Por una parte, los “dos colores” no son simétricos, es decir, con la misma cantidad de material, un array hibridado con uno u otro color sea Cy5 o Cy3, emitirá señales con diferente intensidad. La forma normal de manejar este problema es dye-swapping que consiste en utilizar para una misma comparación dos arrays con dyes cambiados, es decir, si en el primer array la muestra 1 se marca con Cy3 y la muestra 2 con Cy5, con las muestras del segundo array se hace al revés (ver figura 4.3.

Figure 4.3: Representación simplificada de dos diseños

Por otra parte, el hecho de que solo se puedan aplicar dos condiciones a cada array complica el diseño, ya sea porque normalmente hay más de dos condiciones, o porque no es recomendable hibridar directamente dos muestras en un array, creando pares artificiales.

El problema de la asignación eficiente de muestras a microarrays, dado un numero de condiciones a comparar y un número fijo de arrays disponibles, ha sido estudiado de forma exhaustiva.

El diseño utilizado más comúnmente en la comunidad biológica es el diseño de referencia en el que cada condición de interés se compara con muestras tomadas de alguna referencia estándar común a todos las arrays, (ver figura 4.4, imagen (a)).

Figure 4.4: (a) diseño de referencia. (b) diseño en loop.

Los diseños de referencia permiten hacer comparaciones indirectas entre condiciones de interés. La crítica más importante a esta aproximación es que el 50 de las fuentes de hibridación se utilizan para producir la señal del grupo control, de un interés no intrínseco para los biólogos. En contraste, un diseño en loop compara dos condiciones a través de otra cadena de condiciones, por lo que elimina la necesidad de una muestra de referencia.

5 Exploración de los datos, control de calidad y preprocesado

5.1 Introducción

Los estudios realizados con microarrays, sea cual sea la tecnología en que se basan, tienen una característica común: generan grandes cantidades de datos a través de una serie de procesos, que hacen que su significado no siempre sea completamente intuitivo.

Como en todo tipo de análisis, antes de empezar a trabajar con los datos, debemos de asegurarnos de que éstos son fiables y completos y de que se encuentran en la escala apropiada para proporcionar la información que pretendemos obtener de ellos.

En el caso de los microarrays solemos distinguir dos fases previas al análisis de los datos:

Exploración y control de calidad Mediante gráficos y/o resúmenes numéricos estudiamos la estructura de los datos con el fin de decidir si parecen correctos o presentan anomalías que deben ser corregidas.
Preprocesado Incluso si son correctos, los datos “crudos” no sirven para el análisis sino que necesitan preprocesados, lo que puede interpretarse como uno o más de los procesos siguientes:

“Limpiados,” para eliminar la parte de la señal no atribuible a la expresión, llamada “background” o ruído.
“Normalizados” para hacerlos comparables entre muestras y eliminar posibles sesgos técnicos.
Resumidos o “sumarizados” de forma que se tenga un sólo valor por gen.
“Transformados” de forma que la escala sea razonable y facilite el análisis.

A menudo -por un cierto abuso de lenguaje- se denomina “normalización” al preprocesado descrito en las etapas anteriores.

5.1.1 Nivel de análisis y tipo de microarray

Desde la generalización del uso de los microarrays se han desarrollado muchas formas para visualizar los datos y decidir acerca de su calidad. Algunas trabajan sobre los datos obtenidos del escáner, otras lo hacen con los datos normalizados. Algunas sirven tanto para arrays de dos colores como arrays de un color. Otras son específicas de la tecnología. Con el fin de organizar esta multitud de opciones podemos diferenciar:

Datos de bajo nivel Son los datos proporcionados por el escáner contenidos por ejemplo en archivos .gpr (dos colores) o .cel (un color). Estos últimos son binarios por lo que no es posible ni tan sólo visualizarlos sin programas específicos.
Datos de alto nivel Son los datos resultantes del (pre)procesado de los datos de bajo nivel. Básicamente se corresponden con los datos de expresión ya resumidos (“sumarizados”), normalizados o no.

La exploración y el control de calidad pueden basarse en datos de bajo o de alto nivel. En cada caso se pueden aplicar ya sean técnicas generales de visualización de datos, como histogramas, diagramas de caja o visualización en dimensiones reducidas (PCA u otras), o bien técnicas “ad-hoc” para cada tipo de datos como el MA-Plot u otras que se discuten a continuación.

5.1.2 Datos de partida

La estructura de los datos de microarrays de un color o de dos colores difiere considerablemente, tanto a nivel físico (los chips y las imágenes que de ellos se obtiene) como informático, es decir en la forma en que se representan.

5.1.2.1 Arrays de dos colores (o de “cDNA”)

Tradicionalmente los arrays de dos colores o de cDNA se realizaban de forma menos automatizada que los de un color o de Affymetrix. Esto implica que, tras obtener la imagen, el escaneado del archivo “.TIFF” resultante pudiera ser llevado a cabo mediante un software independiente como Genepix, un programa que convierte las imágenes en números y genera un archivo de información (con extensión “.gpr”) a partir del cual pueden calcularse las expresiones relativas, así como valores de calidad para cada “spot” o punto escaneado en la imagen.

Para cada imagen (o sea para cada microarray) hay un archivo “.gpr” que contiene una fila por gen y varias columnas con distintos valores, por ejemplo:

la intensidad para cada canal,
valores resumen de las intensidades y
resultados de controles de calidad (“FLAGS”).

La tabla siguiente muestra lo que serían las primeras filas y columnas de un archivo “.gpr” obtenido mediante el programa Genepix.

Gen	Señal	Fondo	Señal	Fondo	“Flag”	Otros
	Cy5 (R)	Cy5 (Rb)	Cy3 (G)	Cy3 (Gb)
gen-1	3.7547	1.8128	5.0672	1.8496	1	…
gen-2	0.8331	0.9175	1.1536	0.9995	0	…
gen-3	9.8254	0.2781	0.6921	0.5430	1	…
gen-4	9.1539	0.1918	3.8290	0.0014	0	…
gen-5	4.8603	0.2377	0.5338	0.3335	0	…
gen-6	7.8567	1.3941	1.3050	1.4876	1	…
gen-7	2.7619	0.8108	8.4916	1.6518	0	…
gen-8	3.6618	0.1918	0.3770	1.1842	0	…
gen-9	4.4300	0.5888	2.8161	0.8707	1	…
…	…	…	…	…	…	…

Los valores de intensidad se convierten en una única matriz de expresión que contiene una columna por chip con los valores de intensidad relativa obtenidas por ejemplo con una “sumarización” del tipo: \[ \log\frac{R-Rb}{G-Gb}, \] y una fila por gen (mismas filas que archivos “.gpr”).

La tabla siguiente muestra lo que sería una matriz de expresión derivada de cuatro archivos “.gpr” como los de la tabla anterior.

gen	Array-1	Array-2	Array-3	Array-4
gen-1	1.65695	-0.10820	1.69515	8.25137
gen-2	-1.82305	0.11350	1.58807	0.95676
gen-3	0.01561	0.47682	-27.88036	-1.39078
gen-4	0.42709	4.29319	-4.31366	30.96866
gen-5	0.04332	0.24126	-10.27675	0.26680
gen-6	-0.02825	0.68408	2.04163	0.99554
gen-7	3.50549	-8.04635	0.18286	0.22647
gen-8	-0.23261	-1.08477	0.19582	1.11561
gen-9	0.50643	1.52147	0.99092	0.23441
…	…	…	…	…

5.1.2.2 Arrays de un color (Affymetrix)

El resultado de escanear la imagen de un array de affymetrix es un archivo de extensión “.CEL” que, a diferencia de los arrays de dos colores, está en formato binario es decir que solo puede ser leído con programas específicos para ello.

De forma similar a los arrays de dos colores, existe un archivo “.CEL” por cada microarray.

A partir de las intensidades de los archivos “.CEL” se genera la matriz de expresión, que contiene una columna por chip con los valores de intensidad absoluta, y una fila por grupo de sondas. En el caso de arrays de affymetrix existe una gran variedad de algoritmos de “sumarización” y, según cual se utilice, se obtendrá unos u otros valores de expresión. Ahora bien, éstos serán siempre medidas absolutas, es decir independientes del resto de muestras y en una escala arbitraria.

La tabla siguiente muestra las primeras filas de una matriz de expresión sumarizada correspondiente a los primeros genes de uno de los casos resueltos.

ID_REF	GSM188013	GSM188014	GSM188016	GSM188018
1007_s_at	15630.2	17048.8	13667.5	15138.8
1053_at	3614.4	3563.22	2604.65	1945.71
117_at	1032.67	1164.15	510.692	5061.2
121_at	5917.8	6826.67	4562.44	5870.13
1255_g_at	224.525	395.025	207.087	164.835
…	…	…	…	…

5.1.2.3 Datos de ejemplo

Los ejemplos de este capítulo se basarán en dos de los conjuntos de datos descritos en el capítulo 3.

Por un lado, trabajaremos con los datos del conjunto referidos al efecto del tratamiento con estrógenos en cancer de mama descrito en @ref{estrogen}. La exploración se basará en los datos normalizados pero también en los datos “crudos,” de bajo nivel, contenidos en los archivos “.cel.” Estos datos pueden obtenerse directamente del paquete de Bioconductor . Un aspecto interesante de este conjunto de datos es que, junto con 8 muestras “buenas” se proporciona un array defectuoso, denominado “” que facilita el ver como aparecen ciertos gráficos cuando hay problemas.
Por otro lado, usaremos los datos del conjunto , relativos al efecto del tratamiento con CCL4 en la expresión de los hepatocitos. Los datos resultan accesibles al instalar el paquete .

if(!require(BiocManager)) install.packages("BiocManager")
if (!(require(CCl4))){
 BiocManager::install("CCl4")
}
if (!(require(estrogen))){
 BiocManager::install("estrogen")
}
if (!(require(affy))){
 BiocManager::install("affy")
}
if (!(require(affyPLM))){
 BiocManager::install("affyPLM")
}

Aunque, en la actualidad (año 2021) el paquete más utilizado para arrays de este tipo es el paquete oligo este ejemplo utiliza el paquete affy.

library(estrogen)
library(affy)
affyPath <- system.file("extdata", package = "estrogen")
adfAffy = read.AnnotatedDataFrame("phenoData.txt", sep="",  path=affyPath)
affyTargets = pData(adfAffy)
affyTargets$filename = file.path(affyPath, row.names(affyTargets))
affyRaw <- read.affybatch(affyTargets$filename, phenoData=adfAffy)
# show(affyRaw)
actualPath <- getwd()
setwd(affyPath)
allAffyRaw <- ReadAffy()
setwd(actualPath)

El resultado de la lectura es un objeto “affyraw” de clase “affyBatch”:

class(affyRaw)

## [1] "AffyBatch"
## attr(,"package")
## [1] "affy"

print(affyRaw)

## AffyBatch object
## size of arrays=640x640 features (22 kb)
## cdf=HG_U95Av2 (12625 affyids)
## number of samples=8
## number of genes=12625
## annotation=hgu95av2
## notes=

El paquete limma es conocido como paquete para la selección de genes diferencialmente expresados pero también incluye algunas funciones para la lectura y preprocesado de arrays de dos colores.

library("limma")
library("CCl4")
dataPath = system.file("extdata", package="CCl4")
adf = read.AnnotatedDataFrame("samplesInfo.txt", 
    path=dataPath)
#adf
targets = pData(adf)
targets$FileName = row.names(targets)
RG <- read.maimages(targets, path=dataPath, source="genepix")
attach(RG$targets)
newNames <-paste(substr(Cy3,1,3),substr(Cy5,1,3),substr(FileName,10,12), sep="")
colnames(RG$R)<-colnames(RG$G)<-colnames(RG$Rb)<-colnames(RG$Gb)<-rownames(RG$targets)<- newNames
# show(RG)

El resultado de la lectura es un objeto de classe “RGlist.”

class(RG)

## [1] "RGList"
## attr(,"package")
## [1] "limma"

print(RG)

## An object of class "RGList"
## $R
##      DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,]        47        46        40        59        56        40        45
## [2,]      1095        46        39        54        56        40        47
## [3,]        47        46        41        55        52        39        46
## [4,]        46        45        41        56        54        43        47
## [5,]        55        53        45        61        71        47        63
##      DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,]        70        46        65        43        59        46        63
## [2,]        69        45        61        43        62        48        68
## [3,]        81        44        57        41        58        50        63
## [4,]        78        44        60        42        59        48        62
## [5,]        91        51        83        56        71        52        69
##      DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,]        77        40        91        32
## [2,]        71        41        85        32
## [3,]        70        40        39        32
## [4,]        81        41        51        32
## [5,]        81        57        95        37
## 44285 more rows ...
## 
## $G
##      DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,]        77        78        45        45        87        65       128
## [2,]      4558        79        44        47        81        66       131
## [3,]        79        75        46        47        84        66       122
## [4,]        78        82        44        71        86        63       135
## [5,]        93        95        92        57        86        66       132
##      DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,]       107        57        90        99        84        56        68
## [2,]       110        54        89        96       168        54        72
## [3,]       114        54        84        93        78        55        71
## [4,]       113       118        84        97        83        56        75
## [5,]       133        59        95        99       102        64        84
##      DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,]        71        74        83        59
## [2,]        71        75        81        60
## [3,]        73        78        62        59
## [4,]        79        77        70        61
## [5,]        99        83        92        63
## 44285 more rows ...
## 
## $Rb
##      DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,]        41        39        35        47        44        36        40
## [2,]        43        39        35        46        44        36        39
## [3,]        41        40        35        48        44        35        39
## [4,]        40        40        35        49        43        35        39
## [5,]        41        40        36        48        44        35        39
##      DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,]        55        39        49        37        46        42        50
## [2,]        55        37        50        37        46        41        52
## [3,]        59        39        47        37        46        41        52
## [4,]        63        39        48        36        46        41        52
## [5,]        60        39        47        37        46        41        52
##      DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,]        60        36        34        30
## [2,]        59        36        35        31
## [3,]        55        35        35        30
## [4,]        58        36        34        30
## [5,]        56        36        35        30
## 44285 more rows ...
## 
## $Gb
##      DMSCCl319 DMSCCl320 DMSCCl321 DMSCCl329 CClDMS330 CClDMS331 CClDMS332
## [1,]        44        44        32        34        40        39        48
## [2,]        46        43        32        35        41        38        48
## [3,]        43        43        32        35        40        38        49
## [4,]        43        43        32        35        41        38        49
## [5,]        43        42        32        35        41        38        48
##      DMSCCl333 CClDMS379 CClDMS380 CClDMS381 DMSCCl382 DMSCCl384 DMSCCl389
## [1,]        47        35        50        50        45        37        38
## [2,]        48        35        50        49        45        37        38
## [3,]        48        35        51        49        45        37        38
## [4,]        49        35        51        49        46        37        37
## [5,]        49        35        50        49        45        37        38
##      DMSCCl390 CClDMS391 CClDMS393 CClDMS394
## [1,]        38        40        37        36
## [2,]        38        40        37        36
## [3,]        39        40        37        36
## [4,]        38        40        37        36
## [5,]        39        40        37        36
## 44285 more rows ...
## 
## $targets
##            Cy3  Cy5 RIN.Cy3 RIN.Cy5                      FileName
## DMSCCl319 DMSO CCl4     9.0     9.7 251316214319_auto_479-628.gpr
## DMSCCl320 DMSO CCl4     9.0     5.0 251316214320_auto_478-629.gpr
## DMSCCl321 DMSO CCl4     9.0     2.5 251316214321_auto_410-592.gpr
## DMSCCl329 DMSO CCl4     9.0     2.5 251316214329_auto_429-673.gpr
## CClDMS330 CCl4 DMSO     9.7     9.0 251316214330_auto_457-658.gpr
## 13 more rows ...
## 
## $genes
##   Block Row Column           ID            Name
## 1     1   1      1 BrightCorner    BrightCorner
## 2     1   1      2 BrightCorner    BrightCorner
## 3     1   1      3    (-)3xSLv1 NegativeControl
## 4     1   1      4 A_44_P317301        AW523361
## 5     1   1      5 A_44_P386163    NM_001007719
## 44285 more rows ...
## 
## $source
## [1] "genepix"
## 
## $printer
## $ngrid.r
## [1] 1
## 
## $ngrid.c
## [1] 1
## 
## $nspot.r
## [1] 430
## 
## $nspot.c
## [1] 103

5.2 Exploración y control de calidad

Los gráficos son útiles para comprobar la calidad de los datos de microarrays, obtener información sobre cómo se deben preprocesar los datos y comprobar, finalmente, que el preprocesado se haya realizado correctamente.

Siguiendo el esquema presentado en la tabla @ref{tab:tablaDiagnosticos} se presentan a continuación los distintos gráficos utilizados con una breve descripción de lo que representa cada uno y como interpretarlos adecuadamente. El código para generarlos se presenta al final del capítulo como un apéndice.

5.2.1 Control de calidad con gráficos estadísticos generales

5.2.1.1 Histogramas y gráficos de densidad

Estos gráficos permite hacerse una idea de si las distribuciones de los distintos arrays son similares en forma y posición. La figura 5.1 muestra los histogramas correspondientes a los 9 arrays del conjunto de datos .

affySampleNames <- rownames(pData(allAffyRaw))
affyColores <- c(1,2,2,3,3,4,4,8,8)
affyLineas <- c(1,2,2,2,2,3,3,3,3)
hist(allAffyRaw, main="Signal distribution", col=affyColores, lty=affyLineas)
legend (x="topright", legend=affySampleNames , col=affyColores, lty=affyLineas, cex=0.7)

Figure 5.1: La distribución de las expresiones se ve afectada globalmente lo que sugiere más un efecto técnico que biológico

Tanto los gráficos de densidad como los diagramas de cajas de la figura siguiente muestran un desplazamiento global de los valores en todas las muestras del grupo de 48 horas lo que sugiere un efecto técnico, puesto que no es biológicamente plausible que todas los genes se vean más afectados por la mayor tiempo de exposición.

5.2.1.2 Diagramas de caja o “boxplots”

Como los histogramas los diagramas de caja –basados en los distintos cuantiles de las valores– dan una idea de la distribución de las intensidades. La figura 5.2 muestra los diagramas de caja correspondientes a los 9 arrays del conjunto de datos .

boxplot(allAffyRaw, main="Distribución de las expresiones", col=affyColores, las=2)

Figure 5.2: Las distribuciónes de las expresiones se encuentran desplazadas en todas las muestras del grupo de 48

5.2.1.3 Gráficos de componentes principales

El análisis de componentes principales puede servir para detectar si las muestras se agrupan de forma “natural” es decir con otras muestras provenientes del mismo grupo o si no hay correspondencia clara entre grupos experimentales y proximidad en este gráfico. Cuando esto sucede no significa necesariamente que haya un problema pero puede ser indicativo de efectos técnicos -como el conocido efecto “batch”- que podría ser necesario corregir.

plotPCA <- function ( X, labels=NULL, colors=NULL, dataDesc="", scale=FALSE)
{
  pcX<-prcomp(t(X), scale=scale) # o prcomp(t(X))
  loads<- round(pcX$sdev^2/sum(pcX$sdev^2)*100,1)
  xlab<-c(paste("PC1",loads[1],"%"))
  ylab<-c(paste("PC2",loads[2],"%"))
  if (is.null(colors)) colors=1
  plot(pcX$x[,1:2],xlab=xlab,ylab=ylab, col=colors, 
       xlim=c(min(pcX$x[,1])-10, max(pcX$x[,1])+10),
       ylim=c(min(pcX$x[,2])-10, max(pcX$x[,2])+10),
       )
  text(pcX$x[,1],pcX$x[,2], labels, pos=3, cex=0.8)
  titulo <- ifelse(dataDesc=="", "Visualización de las dos primeras componentes", dataDesc)
  title(titulo, cex=0.8)
}

if (!file.exists("datos/affyNorm")){
  allAffyNorm<- rma(allAffyRaw)
  affyNorm <- rma(affyRaw)
  save(allAffyNorm, affyNorm, file="datos/affyNorm.Rda")
}else{
  load(file="datos/affyNorm.Rda")
}

Las figuras 5.3 y 5.4 muestran dos diagramas de componentes principales realizados a partir de los datos normalizados del conjunto de datos . El gráfico de la parte superior que incluye el array defectuoso ilustra que la principal fuente de variabilidad es la diferencia de este array con el resto. Cuando se repite el análisis omitiendo esta muestra puede verse como la principal fuente de variación (eje X) se asocia con el tiempo de exposición (alto a la derecha, bajo (10h) a la izquierda, mientras que la segunda fuente de variación se asocia con la exposición a los estrógenos (alto arriba, bajo abajo).

plotPCA(exprs(allAffyNorm), labels=affySampleNames)

Figure 5.3: La principal fuente de variabilidad si se incluye la muestra defectuos, es la diferencia entre ésta y el resto de observaciones

plotPCA(exprs(affyNorm), labels=colnames(exprs(affyNorm)))

Figure 5.4: Sin la muestra defectuos puede asociarse la primera componente a la diferencia entre 48h y 10h y la segunda a la diferencia ‘high’ vs ‘low’

5.2.1.4 Imagen del chip

Otra forma de ver si las muestras se agrupan según los grupos experimentales, o mediante otros criterios es usando un cluster jerárquico que realiza una agrupación básica de las muestras por grado de similaridad según la distancia que se utilice.

Como en el caso de las componentes principales si las muestras se agrupan según las condiciones experimentales es una buena señal pero si no es así puede deberse a la presencia de otra fuente de variación o bien al hecho de que se trata de un gráfico basado en todo los datos y las condiciones experimentales pueden haber afectado un pequeño número de genes.

La figura 5.5 muestra como se agrupan los datos del conjunto en base a un cluster jerárquico. Como en el caso de las componentes principales tras eliminar el array defectuoso las muestras se separan, primero por el tiempo de exposición y luego por niveles de estrógeno suministrado.

clust.euclid.average <- hclust(dist(t(exprs(affyNorm))),method="average")
plot(clust.euclid.average, labels=colnames(exprs(affyNorm)), main="Hierarchical clustering of samples",  hang=-1, cex=0.7)

Un cluster jerárquico sirve para determinar si las muestras se agrupan de forma natural según los grupos experimentales o si lo hacen de otra forma

Figure 5.5: Un cluster jerárquico sirve para determinar si las muestras se agrupan de forma natural según los grupos experimentales o si lo hacen de otra forma

5.2.2 Gráficos de diagnóstico para microarrays de dos colores

El diagnóstico de arrays de dos canales se basa principalmente en la imagen y en diferentes tipos de gráficos.

5.2.2.1 Diagramas de dispersión y “MA-plots”

La normalización discutida en este mismo capítulo es un punto clave en el proceso de análisis de microarrays y se ha dedicado un gran esfuerzo a desarrollar y probar diferentes métodos (Quackenbush (2002),Yang et al. (2002)). Una razón para ello es que hay diferentes artefactos técnicos que deben ser corregidos para poder ser utilizados, y no cualquier método puede funcionar con todos ellos.

En general, los métodos de normalización se basan en el siguiente principio: _La mayor parte de los genes en un array se pueden expresar o no expresar ante cualquier condición, pero se espera que sólo una pequeña cantidad de genes muestre cambios de expresión entre condiciones}.

Esto da una idea de como debería ser un gráfico de intensidades. Por ejemplo, si no hubiese artefactos técnicos, en arrays de dos canales, una gráfica de dispersión de intensidad del rojo frente al verde debería dejar la mayor parte de los puntos alrededor de una diagonal. Cualquier desviación de esta situación debería ser atribuible a razones técnicas, no biológicas, y por tanto, debería ser eliminada. Esto ha conducido a un método de normalización muy popular consistente en estimar la transformación a aplicar, como una función de las intensidades utilizando el método lowess en la representación transformada de la gráfica de dispersión conocida como el gráfico MA–plot.

La figura 5.6 muestra, en su parte izquierda un gráfico de dispersión del canal rojo frente al verde en un array de dos colores. En concreto se trata del primer array del dataset CCl4 presentado en la introducción.

El hecho de que los datos estén centrados alrededor de la diagonal sugiere que puede no ser necesario normalizar los datos.

Una representación muy popular que ayuda a visualizar mejor esta relación es lo que se conoce como MA-plot, que aparece en la parte derecha (b) de la figura 5.6. Geométricamente representa una rotación del gráfico de dispersión, en la que el significado de los nuevos ejes es:

$A=\displaystyle \frac{1}{2}(\log_2 (R*G))$: El logaritmo de la intensidad media de los dos canales,
$M=\log_2 \displaystyle \frac{R}{G}$: El logaritmo de la expresión relativa entre ambos canales (normalmente conocido como “log–ratio”).

R <- RG$R[,"DMSCCl319"]
G <- RG$G[,"DMSCCl319"]
logR <- log(R)
logG <- log(G)
M <- logR-logG
A <- 0.5*(logR+logG)
opt<- par(mfrow=c(1,2))
plot(logR~ logG, main= "Scatterplot")
abline(lm(logR~logG), col="yellow")
plot(M~A, main= "MA-Plot")
abline(h=0, col="yellow")

Figure 5.6: (a) Gráfico de un canal frente al otro (b) MA-plot (intensidad frente log-ratio)

par(opt)

Si en vez de explorar un array mas moderno como el de la figura anterior lo hacemos con uno de los primeros que se utilizaron, el dataset swirl vemos que el resultado es menos simetrico ,o que sugiere la necesidad de normalizacion.

Este dataset se encuentra disponible en el paquete swirl pero también aparece como ejemplo en el paquete marray. En la viñeta de dicho paquete se realiza una exploración completa de este dataset. Aquí nos limitaremos a reproducir los dos gráficos anteriores para ilustrar un caso en que sí que se pone de manifiesto la necesidad de normalizar los datos.

La clase del objeto “swirl” es distinta a la del objeto “RG” por lo que elk acceso es ligeramente distinto

library(marray)
data(swirl)
R <- swirl@maRf[,1]
G <- swirl@maGf[,1]
logR <- log(R)
logG <- log(G)
M <- logR-logG
A <- 0.5*(logR+logG)
opt<- par(mfrow=c(1,2))
plot(logR~ logG, main= "Scatterplot")
abline(lm(logR~logG), col="yellow")
plot(M~A, main= "MA-Plot")
abline(h=0, col="yellow")

Figure 5.7: (a) Gráfico de un canal frente al otro (b) MA-plot (intensidad frente log-ratio). Estos gráficos sugieren la necesidad de normalizar los datos

par(opt)

Aunque el diagrama de dispersión no lo muestra claramente, en el diagrama “MA” sí que se puede ver como los puntos estan claramente por debajo de lo que sería el eje de simetría horizontal

5.2.2.2 Imagen del array

La imagen del chip (véase la figura ??, izquierda) ofrece una visión rápida de la calidad del array, proporcionando información acerca del balance del color, la uniformidad en la hibridación y en los spots, de si el background es mayor del normal y dela existencia de artefactos como el polvo o pequeñas marcas (rasguños).

5.2.2.3 Histogramas de señales y de la relación señal–ruído

Estos gráficos (véase la figura 5.8, derecha) son útiles para detectar posibles anormalidades o un background excesivamente alto .

Figure 5.8: La relación señal ruido sirve para detectar posibles anormalidades o un background excesivamente alto como medida de calidad

5.2.2.4 Boxplots

Un gráfico muy utilizado es el diagrama de cajas o “boxplot” múltiple con una caja por cada chip. Del alineamiento (o falta de él) y la semejanza (o disparidad) entre las cajas, se deduce si hace falta, o no, normalizar entre arrays.

En el caso de arrays de dos colores pueden utilizarse diagrams de cajas “dentro de arrays” (entre distintos sectores del mismo chip) y “entre arrays.”

5.2.3 Gráficos de diagnóstico para microarrays de un color

5.2.3.1 Imagen del chip

Los arrays de affymetrix contienen millones de sondas por lo que no pueden examinarse a simple vista. A pesar de ello hay diversas formas de obtener una imagen que, en caso de presentar irregularidades pueden indicar algún tipo de problemas como burbujas, arañazos, etc. La figura @ref(c06plotAffy} muestra algunas de las imágenes

La imagen del array de Affymetrix sólo es útil para evidenciar grandes problemas como burbujas, arañazos, etc. En este ejemplo los dos arrays de la izquierda se considerarían aceptables y los de la derecha defectuosos.

Figure 5.9: Imágenes de cuatro microarrays de Affymetrix

5.2.3.2 Gráfico “M-A”

En los chips de dos colores el MA–plot se utliza para comparar los dos canales en cada array (rojo y verde). En cambio, en los chips de Affymetrics, en que sólo hay un canal en cada array, la única forma de definir M (el log ratio) es a partir de la comparación entre pares de de valores, ya sea los arrays dos a dos o bien cada array respecto un valor de referencia que puede ser la mediana, punto a punto, de todos los arrays (véase por ejemplo la figura 5.10).

Figure 5.10: En los chips de Affymetrix la única forma de definir M (el log ratio) es comparar entre diferentes arrays

$M=\log_2(I_1) - \log_2(I_2)$: log ratio $A=\displaystyle \frac{1}{2}(\log_2 (I_1)+\log_2(I_2))$: log de intensidades Donde $I_1$ es la intensidad del array de estudio, e $I_2$ es la intensidad media de arrays. Por lo general, se espera que la distribución en el gráfico se concentre a lo largo del eje M = 0.

5.2.3.3 Modelos de bajo nivel (“Probe-Level-Models” o PLM)

Los modelos de bajo nivel (“Probe-Level-Models” o PLM) ajustan a los valores de intensidad –a nivel de sondas, no de valores totalizados de gen– un modelo explicativo. Los valores estimados por este modelo se comparan con los valores reales y se obtienen los errores o “residuos” del ajuste. El análisis de dichos residuos procede de forma similar a lo que se realiza al analizar un modelo de regresión: Si los errores no presentan ningún patrón especial supondremos que el modelo se ajusta relativamente bien. Si, en cambio, observamos desviaciones de esta presunta aleatoriedad querrá decir que el modelo no explica bien las observaciones, lo cual se atribuirá a la existencia de algún problema con los datos.

Con los valores ajustados del modelo se calculan dos medidas:

La expresión relativa en escala logarítmica ” Relative Log Expression” (RLE) es una medida estandarizada de la expresión. No es de gran utilidad pero debería presentar una distribución similar en todos los arrays.
El error no estandarizado y normalizado o “NUSE” es el más informativo ya que representa la distribucián de los residuos a la que hacíamos referencia más arriba. Si un array es problemático la caja correspondiente en el boxplot aparece desplazada hacia arriba o abajo de las demás.

stopifnot(require(affyPLM))
Pset<- fitPLM(allAffyRaw)

La figura 5.11 muestra los gráficos RLE y NUSE para el conjunto de datos estrogen. En ambos gráficos puede verse como el array defectuoso “” queda claramente diferenciado del resto.

opt<- par(mfrow=c(2,1))
RLE(Pset, main = "Relative Log Expression (RLE)", 
    names=rownames(pData(allAffyRaw)), las=2, cex.axis=0.6)
NUSE(Pset, main = "Normalized Unscaled Standard Errors (NUSE)",  
     names=rownames(pData(allAffyRaw)), las=2, cex.axis=0.6)

Figure 5.11: Graficos de diagnóstico calculados a nivel de sondas PLM

par(opt)

5.3 Normalización de arrays de dos colores

La palabra normalización describe las técnicas utilizadas para transformar adecuadamente los datos antes de que sean analizados. El objetivo es corregir diferencias sistemáticas entre muestras, en la misma o entre imágenes, lo que no representa una verdadera variación entre las muestras biológicas.

Estas diferencias sistemáticas pueden deberse, entre otras, a:

Cambios en la tinción que modifican la intensidad del spot.
La ubicación en el array que puede afectar el proceso de lectura.
Un problem en la placa de origen.
La existencia de diferencias en la calidad de la impresión: pueden presentarse variaciones en los pins y el tiempo de impresión
Camio en los parámetros de la digitalización (escaneo).

A veces puede ser difícil detectar estos problemas , aunque existen algunas formas de saber si es necesaria realizar una normalización. Aqui destacamos dos posibilidades:

Realizar una auto-hibridación. Si hibridamos una muestra con ella misma, las intensidades deberían ser las mismas en ambos canales. Cualquier desviación de esta igualdad, significa que hay un sesgo sistemático.
Detectar artefactos espaciales en la imagen o en la tinción de los gráficos de diagnóstico

5.3.1 Normalización global

Este método esta basado en un ajuste global, es decir en modificar todos los valores una cantidad c, estimada de acuerdo a algún criterio. \[\begin{equation} \log_2 R/G \rightarrow \log_2 R/G-c=\log_2 R/(Gk) \end{equation}\] opciones para $k$ o $c= \log_2k$ son

$c$= mediana o media de log ratio para un conjunto concreto de genes o genes control o genes housekeeping.

La intensidad total de la normalización

$k=\sum R_i/\sum G_i$

5.3.2 Normalización dependiente de la intensidad

En este caso se realiza una modificación específica para cada valor. Esta modificación se obtiene como una función de la intensidad total del gen ($c=c(A)$). \[\begin{equation} \log_2 R/G \rightarrow \log_2 R/G-c(A)=\log_2 R/(Gk(A)) \end{equation}\]

Una posible estimación de esta función puede hacerse utilizando la función lowess (LOcally WEighted Scatterplot Smoothing).

5.4 Resumen y normalización de microarrays de Affymetrix

En los arrays de Affymetrix, como en todos los tipos de microarrays, tras escanear la imagen se obtiene una serie de valores de intensidad de cada elemento del chip. En el caso de estos arrays sabemos que cada valor no corresponde a la expresión de un gen:

Hay múltiples valores (sondas o “probes”) por cada gen, que originan un probeset.
Cada grupo de sondas consiste en múltiples sondas

Originalmente, en los primeros arrays de Affymetrix se introdujo un sistema para estimar el “ruído” y diferenciarlo de la “señal.” Éste consistía en que cada sonda era en ralidad un par de sondas:

Una que coincidía exactamente con el fragmento del gen al que corresponde la sonda y que se denominaba “perfect match.”
Otra, que coincidía con la anterior salvo por el valor central que se sustituía por el nucleótido complementario. Esta segunda sonda se denominaba un “mismatch.”

La idea original consistía en utilizar estos “mismatches’ para tener una medida de hibridación no específica pero en las versiones más recientes se han abandonado, por lo que también se ha abandonado la utilización de los métodos”MAS-4” y “MAS-5” que utilizaban estos valores.

El proceso que convierte las señales individuales en valores de expresión normalizados para cada gen consta de tres etapas:

Corrección del ruido de fondo o “background”
Normalización para hacer los valores comparables
“Sumarización”(Resumen) o concentración de los valores de cada grupo de sondas en un único valor de expresión absoluto normalizado para cada gen.

A menudo los tres pasos se denominan genérica -y erróneamente- “normalización.”

A diferencia de los chips de ADNc, aquí las medidas de expresión son absolutas (no se compara una condición contra otra) dado que cada chip se hibrida con un única muestra.

Hay muchos métodos para estimar la expresión (más de 30 publicados). Cada método contempla de forma explícita o implícita las tres formas de preprocesado: corrección del fondo, normalización y resumen.

Aquí nos centraremos únicamente en el método más popular que continua siendo utilizado trece años después de su introducción: El “Robust Multi-chip Average” o RMA

5.4.1 El método RMA (Robust Multi-Array Average)

Para compensar algunas deficiencias de los primeros métodos de resumen y normalización de arrays de Affymetrix, Irizarry y sus colegas introdujeron en 2003 (~) un método basado en la modelización de las intensidades de las sondas que, en vez de basarse en las distintas sondas de un gen dentro de un mismo array se basa en los distintos valores de la misma sonda entre todos los arrays disponibles,

Esquemáticamente los pasos que realiza este método son:

Ajusta el ruído de fondo (background) basándose solo en los valores PM y utilizando un modelo estadístico complejo en el que combina la modelización de la señal mediante una distribución exponencial con la del ruído mediante una distribución normal.
Toma logaritmos base 2 de cada intensidad ajustada por el background.
Realiza una normalización por cuantiles de los valores del paso 2 consistente en substituir cada valor individual por el que tendría la misma posición en la distribución empírica estimada sobre todas las muestras, es decir los promedios de las distribuciones de los valores ordenados de cada array (véase figura )
Estima las intensidades de cada gen separadamente para cada conjunto de sondas. Para ello realiza una técnica similar a una regresión robusta denominada “pulido de medianas” (median polish) sobre una matriz de datos que tiene los arrays en filas y los grupos de sondas en columnas.

Como resultado final de todos los pasos anteriores se obtiene la matriz con los datos sumarizados y normalizados. A pesar de no estar exento de críticas como la que afirma que este procedimiento “compacta” los valores reduciendo su variabilidad natural, este método se ha convertido en el estándar “de facto” actualmente por muchos usuarios de Bioconductor.

Figure 5.12: El método RMA incluye una normalización por cuantiles como la representada esquemáticamente en esta figura

5.5 Filtraje no específico

El objetivo del filtraje es eliminar aquellos genes o aquellos spots cuyas imágenes o señales sean erróneas por diferentes motivos, con el fin de reducir el ruído de fondo.

Los principales procesos de filtrado son:

Eliminación de los spots marcados como erróneos mediante flags y que son debidos a problemas en la hibridación o en el escaneo.
Eliminación de spots o genes con señales muy bajas debido a problemas en el spotting o a que no ha habido hibridación en ese spot (Filtraje por señal).
Eliminación de genes que no presenten una variación significativa en su señal, entre distintas condiciones experimentales (Filtraje por variabilidad).

Ante la duda, se suele aconsejar ser conservador y reducir la operación de filtraje al mínimo, puesto que cuesta mucho convencerse, y convencer a los investigadores, de que el beneficio del filtraje compensa la posibilidad de eliminar, por error, un gen o una señal interesantes.

6 Selección de genes diferencialmente expresados

6.1 Introducción

El motivo más habitual para el que se suelen utilizar microarrays es la búsqueda de genes cuya expresión cambia entre dos o más condiciones experimentales, por ejemplo a consecuencia de un tratamiento, una enfermedad u otras causas (distintos tiempos, distintas lineas celulares, …).

El problema consiste en identificar estos genes y suele denominarse selección de genes diferencialmente expresados (“DEG”) o bien comparación de clases.

El problema de seleccionar genes diferencialmente expresados se traduce de manera casi inmediata al problema estadístico de comparar variables y, en años recientes, se han desarrollado un gran número de métodos estadísticos para resolverlo. La mayoría son extensiones de los métodos estadísticos clásicos, pruebas-t o análisis de la varianza, adaptados en uno u otro sentido para tener en cuenta las peculiaridades de los microarrays.

Aunque el problema de la selección de genes diferencialmente expresados puede relacionarse directamente con la realización de pruebas estadísticas, en el caso de los microarrays, el hecho de que haya dos tecnologías que miden la expresión de dos formas distintas hace que se deba diferenciar la metodología a emplear en cada caso. Los arrays de dos colores combinan dos muestras en un chip y generan una medida de expresión relativa. Esto hace que para comparar dos muestras de un mismo individuo sean la opción naturalmente más apropiada

En el caso de querer comparar muestras independientes de diferentes individuos los arrays de un color son la mejor opción. Evidentemente, lo más común será disponer de una sola técnica y tener que adaptar los análiaisis estadísticos a la misma.

Vamos a plantear un posible esquema de trabajo, para situar la mejor opción en cada caso:

Situación 1: experimento con 5 individuos diabéticos y 5 no diabéticos, independientes entre si (muestras independientes)
- Caso 1: Arrays de cDNA (2 colores): Utilizaríamos 5 arrays Diabético/Referencia y 5 arrays No diabético/Referencia
- Caso 2: Arays de Affymetrix (1 color): Utilizaríamos 5 arrays de Diabético y 5 arrays de No diabético
Situación 2: experimento con 6 individuos de los que se ha tomado una muestra de tejido sano y otra de tejido tumoral (muestras apareadas o dependientes)
- Caso 3: Arrays de cDNA (2 colores): Utilizaríamos 6 arrays, uno por individuo, y en cada uno se realizaría la comparación Tejido Tumoral/Tejido sano.
- Caso 4: Arays de Affymetrix (1 color): Utilizaríamos 12 arrays, 2 por individuo, 6 con muestras de Tejido Tumoral y 6 con muestras de Tejido sano.

6.2 Selección de genes diferencialmente expresados

6.2.1 Medidas naturales para comparar dos muestras

Recordemos que una vez se han hecho los experimentos con microarrays y obtenida la señal, los datos que se disponen son los logaritmos del valor detectado por el escaner.Esto hace que algunas operaciones que se realicen tengan en cuenta esta característica. Segun si la comparación a realizar se llevará a cabo con datos independientes (2 muestras, casos 1 y 2) o con datos dependientes (muestras apareadas, casos 3, 4) algunas medidas naturales o razonables para la comparación de expresiones son las siguientes:

Para comparaciones directas, con expresiones relativas entre muestras apareadas o bien diferencias apareadas de expresiones absolutas:
- log ratio promedio: $\overline{R}=\frac 1n \sum_{i=1} R_i$
- t-test de una muestra $\frac{\overline{R}}{SE}$, donde SE estima el error estándar del log ratio promedio
- t-test robusto: Substituir en el anterior medidas robustas del error estándar
Para comparaciones indirectas entre muestras independientes de expresiones relativas o absolutas:
- Diferencia media $\overline{R}_1-\overline{R}_2= \frac 1n_1 \sum_{i=1} ^{n_1}R_i- \frac 1n_2 \sum_{j=1} ^{n_2}R_j$
- t-test (clásico) de dos muestras $\frac{\overline{R}_1-\overline{R}_2} {SE_{12}\sqrt{\frac 1n_1 +\frac 1n_2 }}$
- t-test robusto de dos muestras: Substituir en el anterior medidas robustas del error estándar

6.2.2 Un primer ejemplo

Consideremos la tabla siguiente que representa una matriz de expresión simplificada que contiene las expresiones relativas (por ejemplo entre tejido tumoral y sano del mismo individuo) de 5 genes en 6 muestras.

Podemos calcular las medidas descritas para el caso de una muestra para decidir si un gen está expresado o no lo está. Se discutirá más adelante como precisar esto pero de momento nos quedaremos con la idea de que si la medida escogida es (cercana a) cero el gen no está diferencialmente expresado y si es mayor o menor que cero si que lo está. Nos referimos a “cero” porque estamos hablando de logaritmos de razones: si la expresión es la misma en ambas condiciones el cociente es uno y su logaritmo es cero.

Vale la pena insistir en el concepto de expresión diferencial: no nos preocupa cual es la expresión del gen en una u otra muestra sinó si son distintas.

La tabla anterior sugiere que podría considerarse el gen A está diferencialmente expresado (promedio y t–test altos) mientras que el gen B o el D no lo están (promedio y test-t próximos a cero). Los genes C y D pueden llevar a conclusiones contradictoria según nos basemos en el promedio o el test t.

Si se observan los valores del test t del gen C se concluye que el gen no aparenta estar diferencialmente expresado. Si en cambio se observa su promedio parece que si que lo esté.

En el gen E pasa exctamente lo opuesto. La explicación de estas aparentes contradicciones se halla en el error estándar. En el gen C es muy elevado, debido a que el valor (20) es probablemente un “outlier.” En el gen D el error estándar es muy bajo por lo que, al encontrarse en el denominador del t-test aumenta artificialmente su valor.

6.2.3 Selección de genes usando tests estadísticos

Vamos a plantear la forma como se aborda este problema en el estudio de datos de microarrays. Dadas las características propias de este tipo de datos se consideran otras formas de estimar el error estándar que no sean tan sensibles a valores extremos o muy bajos.

Consideremos el gen $g$. Si llamamos:

$R_g$ log-ratio medio observado.
$SE_g$ error estándar de $R_g$ estimado a partir de los datos en el gen $g$ .
$SE$ error estándar de $R_g$ estimado a partir de los datos con la información de todos los genes.

Podemos considerar dos variantes para el test-t:

Test-t global: se calcula en base a un único estimador de SE para todos los genes: \[t=R_g/SE,\]
Test-t específico: Utiliza un estimador distinto del error estandar para cada gen: \[t=R_g/SE_g.\] Cada aproximación tiene sus pros y sus contras como muestra la tabla siguiente: \begin{center}

Test

Test-t Global

Test-t específico

\end{center}

En la práctica muchos métodos de selección de genes diferencialmente expresados han acabado buscando un compromiso entre ambas aproximaciones para lo que proponen o derivan fórmulas que de alguna forma ponderan o combinan dos estimaciones del error estándar: una basada en todos los genes y otra específica de cada gen. La tabla siguiente ilustra como algunos de los métodos más utilizados en la bibliografía incorporan esta idea.

Al hecho de incluir un coeficiente que tenga en cuenta la variabilidad de todos los genes en el array para estimar el error estándar de cada gen se le denomina moderación de la varianza (“variance shrinkage”) y es una de las aproximaciones en que existe cierto consenso (Allison et al. (2006)) acerca de que sirven para mejorar la selección de genes diferencialmente expresados.

6.2.3.1 Ejemplo de utilización del test-t

Como ejemplo utilizaremos el conjunto de datos celltypes y supondremos que disponemos ya de los datos normalizados y filtrados almacenados en un objeto expressionSet.

Este proceso es sencillo de hacer siguiendo los pasos del capiítulo anterior. Para simplificar el ejemplo una vez cargados y normalizados filtraremos los datos usando funcion nSFilter de la librería genefilter con sus opciones por defecto, lo que nos dejara un total de 10254 sondas en vez de las 45101 originales.

library(affy)
library(genefilter)
affyPath<- "datos/celltypes/celfiles"
cellTypesTargets<-  read.AnnotatedDataFrame("targets.txt", sep="\t",  path=affyPath)
cellTypesTargets$filename = file.path(affyPath, row.names(cellTypesTargets))
cellTypesRaw <- read.affybatch(cellTypesTargets$filename, phenoData=cellTypesTargets)
eset_rma <- rma(cellTypesRaw)
Filtered <- nsFilter(eset_rma)
eset_rma_filtered <- Filtered[["eset"]]
save(eset_rma, eset_rma_filtered, file="datos/celltypes/celltypes-normalized.rma.Rda")

Si consideramos que el campo treat del objeto contiene la información de los grupos a comparar (ratones estimulados con LPS frente a los que no han sido tratados), podemos realizar un primer análisis utilizando el test-t. Para ello utilizaremos el la función rowttests del paquete genefilter que realiza un test $t$ sobre cada una de las filas (genes) de una matriz de expresión.

stopifnot(require(Biobase))
load (file="./datos/celltypes/celltypes-normalized.rma.Rda")
my.eset <- eset_rma_filtered
grupo_1 <- as.factor(pData(my.eset)$treat)
stopifnot(require(genefilter))
teststat <-rowttests(my.eset, "treat")
print(teststat[1:5,])

##                statistic         dm      p.value
## 1424378_at   -10.6739944 -1.1113660 8.715342e-07
## 1417675_a_at  20.7390631  0.9451576 1.504712e-09
## 1436530_at     9.2313983  1.4943922 3.291690e-06
## 1428603_at    -3.4882096 -0.4013432 5.840466e-03
## 1458888_at     0.7083721  0.2343463 4.948949e-01

Cuanto mayor sea el valor absoluto del estadístico $t$ mayor es la probabilidad de que el gen esté diferencialmente expresado.

6.3 Significación estadística y expresión diferencial

Como hemos visto en la sección anterior dos medidas naturales para la selección de genes son el promedio de “log-ratios,” o la diferencia de promedios en el caso de muestras independientes, o el valor del estadístico de test ($t$-test) de una o dos muestras según si se trata de muestras apareadas o independientes respectivamente. Los primeros estudios de microarrays eran muy costosos y se hacían con pocas o incluso ninguna réplica por condición experimental. En estas situaciones la única forma fiable de detectar una diferencia de expresión era a traves del “log-ratio” o su diferencias.

Rápidamente se puso en evidencia que para poder obtener los genes que estaban realmente diferencialmente expresados era preciso disponer de un soporte estadístico que permitiera tener en cuenta la variación aleatoria existente entre muestras.

En la práctica esto se reduce a afirmar que si, además de la diferencia de expresión entre las condiciones experimentales, se lleva a cabo un test estadístico dispondremos de una medida objetiva, el p–valor que nos servirá para decidir qué genes se declaran diferencialmente expresados, a saber, aquellos en los que el p–valor del test sea inferior a un cierto umbral como 0.05 o 0.01.

Como es sabido, un test estadístico procede decidiendo rechazar la hipótesis nula si el p–valor es más pequeño que el nivel de significación del test.

Siguiendo con el ejemplo anterior podemos ordenar los resultados de los tests en base a los p–valores:

ranked <-teststat[order(teststat$p.value),]
print(ranked[1:5,])

##              statistic        dm      p.value
## 1449383_at   -48.61014 -2.044928 3.276616e-13
## 1451421_a_at -40.72173 -1.445182 1.909283e-12
## 1415929_at   -39.77036 -0.812879 2.415238e-12
## 1450826_a_at  37.62239  5.147992 4.193941e-12
## 1448303_at   -31.92365 -2.283765 2.140612e-11

Ahora podríamos seleccionar, por ejemplo los genes cuyo p–valor fuera inferior a 0.01

selectedTeststat <- ranked[ranked$p.value < 0.01,]

Esto deja un total de {r} dim(selectedTeststat)[1] con un p–valor inferior a 0.01

6.3.1 “Volcano plots”

Si se opta por computar los valores de significación (p–valores) de los genes, resulta interesante comparar el tamaño del cambio del nivel de significación estadístico. El “volcano plot” es una representación gráfica que permite ordenar los genes a lo largo de dos dimensiones, la biológica, representada por el “fold change” y la estadística representada por el logaritmo negativo del p–valor.

En la escala horizontal se representa el cambio entre los dos grupos (en escala logarítmica, de manera que la regulación positiva o negativa se representa de forma simétrica). En la escala vertical se representa el p–valor del test en una escala logarítmica negativa, de forma que los p–valores más pequeños aparecen mayores.

Así pues puede considerarse que el primer eje indica impacto biológico del cambio (a más efecto biológico mayor “fold-change”) y el segundo muestra la evidencia estadística, o la fiabilidad del cambio.

Como veremos más adelante en esta sección, el paquete limma permite realizar volcano-plots de forma sencilla a partir de los resultados de un anàlisis, pero para hacer un volcano plot tan sólo se necesita una lista de p-valores y los efectos (“fold-change”) asociados.

La figura (ref?)(fig:volcano0) muestra un “volcano-plot” para este ejemplo de los “celltypes” para la comparación “LPS” frente a “Medium.”

FC <- teststat$statistic
pVal<-teststat$p.value
Y <- -log(pVal)
plot(Y~FC)

También podemos adaptar alguna función “ad-hoc” como la siguiente, tomada de :

library(ggrepel)
# plot adding up all layers we have seen so far
volcanoP<- function (de,log2FoldChange,  pvalue, 
                     diffexpressed=NULL, col=NULL, delabel=NULL){
  ggplot(data=de, aes(x=log2FoldChange, y=-log10(pvalue), col=diffexpressed, label=delabel)) +
        geom_point() + 
        theme_minimal() +
        geom_text_repel() +
        scale_color_manual(values=c("blue", "black", "red")) +
        geom_vline(xintercept=c(-0.6, 0.6), col="red") +
        geom_hline(yintercept=-log10(0.05), col="red")
}
diffexpressed <- ifelse(abs(teststat$statistic)>10, TRUE, FALSE)
label <- rep(NA, length(diffexpressed))
label[diffexpressed] <- rownames(teststat)[diffexpressed]
volcanoP (de=teststat, log2FoldChange=teststat$statistic,  pvalue=teststat$p.value,
          diffexpressed = diffexpressed, delabel = label)

6.3.2 Potencia y tamaño muestral

Tal como se ha indicado más arriba, para realizar un buen test suele controlarse la probabilidad de error de tipo I (de falsos positivos) y buscar, de entre los tests candidatos, aquellos que tengan una menor probabilidad de error de tipo II, o equivalentemente una mayor potencia. A partir de este planteamiento existe, en el contexto estadístico estándar, multitud de formas de determinar cual debe ser el tamaño muestral necesario para obtener una potencia dada fijados el tamaño de efecto (“fold-change”) y la probabilidad de error de tipo I.

En el caso de los microarrays se han desarrollado diversas fórmulas para realizar cálculos de este tipo, pero la compleja estructura de los datos microarrays hace que sean relativamente discutibles.

Simon (Simon et al. (2003)) sugiere la fórmula siguiente que es una generalización de las fórmulas clásicas para problemas de dos muestras:

El tamaño total requerido para detectar genes diferencialmente expresados en al menos una diferencia $\delta$ con una probabilidad de error de tipo I (FP), $\alpha$ y una probabilidad de error de tipo II (FN) $1-\beta$ se calcula: \[ n=\frac{4(z_{\alpha/2}+z_\beta)^2}{(\delta/\sigma)^2}, \] donde $z_\alpha$ y $z_\beta$ son los percentiles 100$\alpha/2$ y 100$\beta$ de la distribución Normal $N(0,1)$ y $\sigma$ es la desviación estandar de un gen dentro de una clase (de un grupo). Obviamente $\sigma$ es siempre desconocida por lo que, sin una muestra piloto con que estimarla el calculo e más imaginativo que realista.

Además de esto, el número de arrays usualmente recomendado queda lejos de la cantidad asequible para la mayor parte de los experimentos (Lee and Whitmore (2002),Tibshirani (2006)). Lo que muchos usuarios hacen a la práctica, es buscar un equilibrio entre los costes y la reproducibilidad y, de hecho, tienden a usar una cantidad fija de arrays tal como 3 o 5 sin consideraciones adicionales.

Por ejemplo si ponemos $\alpha=0.001$, $1-\beta=0.95$, $\delta=1$ y estimamos $\sigma$ entre todas las muestras el número de réplicas biológicas que necesitaremos será de 35.8.

6.4 El problema de la múltiplicidad de tests (“multiple testing”)

El análisis de microarrays se realiza en base gen a gen e involucra múltiples tests, miles probablemente. Esto significa que, a medida que crece el número de genes, la probabilidad de declarar erroneamente al menos un gen diferencialmente expresado va en aumento, y si no se realiza algún tipo de ajuste el número de falsos positivos será tanto más alto cuantos más genes estemos analizando.

Hay muchas formas de intentar controlar estas probabilidades de error y puede verse un excelente revisión en Dudoit Dudoit, Shaffer, and Boldrick (2003)).

De forma simplificada consideramos las dos aproximaciones más utilizadas.

Una posibilidad es mirar de controlar la probabilidad de obtener al menos un falso positivo o “Family-wise-error-rate (FWER).” El más popular de estos métodos de control es la corrección de Bonferroni, consistente en multiplicar el p–valor por el número de tests realizados. La misma Dudoit y muchos otros autores han desarrollado variantes de los métodos clásicos de ajuste FWER usando por ejemplo tests de permutaciones.

El criterio FWER es quizás demasiado restrictivo dado que el control de los falsos positivos implica un considerable incremento de falsos negativos. En la práctica, sin embargo, muchos biólogos parecen estar dispuestos a aceptar que se produzcan algunos errores, siempre y cuando esto permita realizar descubrimientos. Por ejemplo un investigador debe considerar aceptable cierta pequeña proporción de errores (digamos del 10 al 20%) entre sus descubrimientos. En este caso, el investigador está expresando interés en controlar el porcentaje de falsos descubrimientos (FDR), es decir lo que es la proporción de falsos positivos sobre el total de genes inicialmente identificados como expresados diferencialmente. A diferencia del nivel de significación que queda determinado antes de examinar los datos, la FDR es una medida de confianza a posteriori ya que emplea información disponible en los datos para estimar las proporciones de falsos positivos que han tenido lugar. Si se obtiene una lista de los genes expresados diferencialmente en los que el FDR se controla hasta, digamos, el 20%, cabe esperar que el 20% de estos genes representen falsos positivos. Lo cual supone un enfoque menos restrictivo que controlar el FWER.

La decisión de controlar el FDR o el FWER depende de los objetivos del experimento. Si, por ejemplo, el objetivo es la captura de genes es razonable permitir cierta cantidad de falsos positivos y es preferible seleccionar FDR. Si por el contrario se trabaja con una lista de un tamaño menor al deseado para verificar la expresión de ciertos genes específicos, entonces el FWER es el criterio apropiado.

El ejemplo siguiente muestra como se realiza el ajuste de p–valores usando el paquete multtest de Bioconductor para ajustar por los métodos de Bonferroni (FWER), Benjamini & Hochberg (FDR) o Benjamini & Yekutieli (BY, FDR).

stopifnot(require(multtest))
procs <- c("Bonferroni","BH", "BY")
adjPvalues <- mt.rawp2adjp(teststat$p.value, procs)
names(adjPvalues)

## [1] "adjp"    "index"   "h0.ABH"  "h0.TSBH"

ranked.adjusted<-cbind(ranked, adjPvalues$adjp)
head(ranked.adjusted)

##              statistic        dm      p.value         rawp   Bonferroni
## 1449383_at   -48.61014 -2.044928 3.276616e-13 3.276616e-13 3.359842e-09
## 1451421_a_at -40.72173 -1.445182 1.909283e-12 1.909283e-12 1.957778e-08
## 1415929_at   -39.77036 -0.812879 2.415238e-12 2.415238e-12 2.476585e-08
## 1450826_a_at  37.62239  5.147992 4.193941e-12 4.193941e-12 4.300467e-08
## 1448303_at   -31.92365 -2.283765 2.140612e-11 2.140612e-11 2.194983e-07
## 1418645_at   -28.74611 -4.126867 6.044555e-11 6.044555e-11 6.198087e-07
##                        BH           BY
## 1449383_at   3.359842e-09 3.296908e-08
## 1451421_a_at 8.255283e-09 8.100651e-08
## 1415929_at   8.255283e-09 8.100651e-08
## 1450826_a_at 1.075117e-08 1.054978e-07
## 1448303_at   4.389966e-08 4.307737e-07
## 1418645_at   1.033014e-07 1.013665e-06

Si seleccionamos los genes en base a su p–valor ajustado por ejemplo por le método de Banjamini y Yekutieli se obtienen los siguientes genes

selectedAdjusted<-ranked.adjusted[ranked.adjusted$BY<0.001,]
nrow(selectedAdjusted)

## [1] 420

6.5 Modelos lineales para la selección de genes: `limma`

En la sección anterior se ha discutido el uso del test $t$ y sus extensiones para la selección de genes diferencialmente expresados en situaciones relativamente sencillas, es decir cuando sólo hay dos grupos.

En muchos estudios el número de grupos a considerar es más de dos y las fuentes de variabilidad pueden ser más de una, por ejemplo una puede ser el tratamiento pero otra puede ser la edad de los individuos o cualquier otra condición fijada por el investigador o derivadad de la heterogeneidad de las muestras.

En estas situaciones una aproximación razonable en problemas con una variable respuesta es el análisis de la varianza, discutido en el capítulo here. En esta sección se presenta una aproximación equivalente que de forma general se denomina el modelo lineal. Este método -que engloba el análisis d la varianza y la regresión- es uno de los más usados en estadística y ha sido popularizado en el campo de microarrays gracias a los trabajos de Gordon Smyth quien ha creado el paquete limma que se ha convertido en la herramiente más utilizada para el análisis de microarrays.

6.5.1 El modelo lineal general

El modelo lineal (ver por ejemplo Faraway, Faraway (2004)) es un marco general para la modelización y el análisis de datos estadística.

EL método consiste en asumir una relación lineal entre los valores observados de una variable respuesta y las condiciones experimentales. A partir de aquí se obtienen estimadores para los parámetros del modelo y de sus errores estándar, y (con algunas condiciones extra) es posible hacer inferencia acerca del experimento.

La aplicación de modelos lineales puede ser visto como un proceso secuencial, con los siguientes pasos:

Especificar el diseño del experimento: qué muestras se asignan a qué condiciones.
(Re-)Escribir un modelo lineal para este diseño en forma de $\mathbf{Y=X\beta+\varepsilon}$, donde $\mathbf{X}$ se denomina la matriz de diseño.
Una vez que el modelo se ha especificado aplicar la teoría general de estimar los parámetros y los contrastes (comparaciones entre los valores de los parámetros).
Si se cumplen ciertas condiciones de validez para el modelo es posible realizar inferencia sobre los parámetros del modelo, es decir se pueden contrastar hipótesis sobre dichos parámetros.

El esquema anterior se puede aplicar a casi cualquier tipo de situación experimental. En la sección siguiente se presentan un par de ellas.

6.5.2 Ejemplos de situaciones modelizables linealmente

6.5.2.1 Ejemplo 1: Experimento “Swirl–Zebrafish”

Swirl es una mutación puntual que provoca defectos en la organización del embrión en desarrollo a lo largo de su eje dorsal-ventral. Como resultado, algunos tipos celulares se reducen y otros se expanden. Un objetivo de este trabajo fue identificar los genes con expresión alterada en el mutante Swirl en comparación con “wild zebrafish.”

El diseño experimental para este estudio fue el siguiente:

Array	Cy3	Cy5
1	W	M
2	M	W
3	W	M
4	M	W

Cada microarray contenía 8848 sondas de cDNA (genes o sequencias EST).
Cuatro réplicas por array (slide): 2 juegos de pares de intercambio de color
El cDNA del mutante swirl ($S$) se marca con Cy5 o Cy3 y el cDNA del “wild type” se marca con el otro

El modelo lineal derivado del diseño anterior fue:

El parámetro de interés es: $\alpha= \mathbf{E}\left (log\frac{S}{W}\right )$ .
Las muestras 1 y 3 estan marcadas con : S (Verde=“Green”) y W (Rojo=“Red”), y las muestras 2 y 4 son “dye-swapped.”
El modelo, $\mathbf{y=X\alpha+\varepsilon}$, es:

\[ \begin{array}{ccccc} y_1 &=&\alpha&+&\varepsilon _1 \\ y_2 &=&-\alpha&+&\varepsilon _2 \\ y_3 &=&\alpha&+&\varepsilon _3 \\ y_4 &=&-\alpha&+&\varepsilon _4\\ \end{array} \Longrightarrow \left( \begin{array}{r} y_1 \\ y_2 \\ y_3\\ y_4 \end{array} \right) = \underbrace{\left( \begin{array}{r} 1 \\ -1 \\ 1\\ -1 \end{array} \right)}_{\text{Matriz de Diseño}, \mathbf{X}} \alpha + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \end{array} \right) \]

6.5.3 Ejemplo 2: Comparación de tres grupos

Los plásmidos IncHI codifican genes de resistencia múltiple a los antibióticos en S. enterica.

El plásmido R27 de la cepa salvaje es termosensible al transferirse.

Algunos fenotipos mutantes relacionados con hha y hns cromosómicos participan en diferentes procesos metabólicos de interés en la conjugación termoregulada.

El objetivo del experimento es encontrar qué genes se expresean diferencialmente en tres tipos de mutantes diferentes, $M_1$, $M_2$ y $M_3$.

Este experimento debe ser planteado de forma diferente según el tipo de arrays (uno o dos colores) y qué comparaciones son las de mayor interés.

Array de dos colores
- A diseño de referencia: Hibridar cada Mutante ($M_i$) vs. Salvaje (“Wild type”) ($W$).
- A diseño loop: Hibridar cada mutante el uno al otro en un doble “loop” que incluye (“dye-swapping”).
Array de un color: hibridar mutantes y “wild types” separadamente.

Permite la comparación directa de Mutantes vs “Wild.”
Número de parámetros a estimar es 3, relación intuitiva entre el número de parámetros y mutantes.
Las comparaciones de Mutante vs Mutante son menos eficientes.

Permite la comparación directa de Mutantes vs Mutantes.
El número de parámetros a estimar es 2. Menos intuitivo.
Las comparaciones Mutante vs Mutante son más eficientes.

Permite la comparación directa de
- Mutante vs “Wild”
- Mutante vs Mutante
El número de parámetros a estimar es 4.
Todas las comparaciones se pueden hacer de forma eficiente.

Modelo, $\mathbf{y=X\alpha+\varepsilon}$, y contrastes $\mathbf{C'\beta}$

Parámetros del modelo: \[\alpha_1= \mathbf{E}\left (log\frac{M_1}{W}\right ),\ \alpha_2= \mathbf{E}\left (log\frac{M_2}{W}\right ),\ \alpha_3=\mathbf{E}\left (log\frac{M_3}{W}\right ).\]
Contrastes: Comparaciones interesantes. \[\begin{array}{ccc} \beta_1 &=& \alpha_1-\alpha_2,\\ \beta_2 &=& \alpha_1-\alpha_3,\\ \beta_3 &=& \alpha_2-\alpha_3. \end{array}\]

\[ \left( \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \end{array} \right) = \underbrace{\left( \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ -1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & -1 \end{array} \right)}_{\text{Matriz de Diseño}, \mathbf{X}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \end{array} \right) + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6 \end{array} \right) \] \[ \left( \begin{array}{r} \beta_1 \\ \beta_2 \\ \beta_3 \end{array} \right) = \underbrace{\left( \begin{array}{rrr} 1 & -1 & 0 \\ 1 & 0 & -1 \\ 0 & 1 & -1 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \end{array} \right). \]

Modelo, $\mathbf{y=X\alpha+\varepsilon}$, y contrastes $\mathbf{C'\beta}$

Parámetros del modelo: \[ \alpha_1= \mathbf{E}\left (log\frac{M_1}{M_2}\right ),\ \alpha_2= \mathbf{E}\left (log\frac{M_2}{M_3}\right ). \] $\alpha_3$ no es necesaria: $\log \left( \frac{M_1}{M_3}\right )=\log\left(\frac{M_1}{M_2}\right )-\log\left(\frac{M_2}{M_3}\right)$.
Contrastes: Algunas de las comparaciones deseadas son precisamente los parámetros. \[\begin{array}{ccc} \beta_1 &=& \alpha_1,\\ \beta_2 &=& \alpha_2,\\ \beta_3 & =& \alpha_1+\alpha_2. \end{array}\]

\[ \left( \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \end{array} \right) = \underbrace{\left( \begin{array}{rr} 1 & 0 \\ 0 & 1 \\ 1 & -1 \\ -1 & 0 \\ 0 & -1 \\ -1 & 1 \end{array} \right)}_{\text{Matriz de Diseño}, \mathbf{X}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \end{array} \right) + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6 \end{array} \right) \]

\[ \left( \begin{array}{r} \beta_1 \\ \beta_2 \\ \beta_3 \end{array} \right) = \underbrace{\left( \begin{array}{rrr} 1 & 0 \\ 0 & 1 \\ 1 & +1 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \end{array} \right). \]

Modelo, $\mathbf{y=X\alpha+\varepsilon}$, y contrastes $\mathbf{C^{1'}\beta}$, $\mathbf{C^{2'}\beta}$

Parámetros del modelo: \[\alpha_1= \mathbf{E} (log{M_1} ),\ \alpha_2= \mathbf{E} (log{M_2} ),\ \alpha_3=\mathbf{E} (log{M_3} ), \alpha_4=\mathbf{E} (log{W} ). \]
Contrastes: Dos posibles conjuntos de comparaciones interesantes.

Comparación entre tipo de mutantes $(\mathbf{C^{1'}\beta})$ \[\begin{array}{ccc} \beta_1^1 &=& \alpha_1-\alpha_2,\\ \beta_2^1 &=& \alpha_3-\alpha_2,\\ \beta_3^1 &=& \alpha_2-\alpha_3. \end{array}\]
Comparación entre cada mutantes y el wild type $(\mathbf{C^{2'}\beta})$ \[\begin{array}{ccc} \beta_1^2 &=& \alpha_4-\alpha_1,\\ \beta_2^2 &=& \alpha_3-\alpha_1,\\ \beta_3^2 &=& \alpha_2-\alpha_1. \end{array}\]

\[ \left( \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \\ y_7 \\ y_8 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \end{array} \right)}_{\text{Matriz de Diseño}, \mathbf{X}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right) + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6\\ \varepsilon_7\\ \varepsilon_8 \end{array} \right) \]

\[ \left( \begin{array}{r} \beta_1^1 \\ \beta_2^1 \\ \beta_3^1 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} 1 & -1 & 0 & 0\\ 1 & 0 & -1 & 0\\ 0 & 1 & -1& 0 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C^1}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right). \]

\[ \left( \begin{array}{r} \beta_1^2 \\ \beta_2^2 \\ \beta_3^2 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} 1 & 0 & 0 & -1\\ 0 & 1 & 0 & -1\\ 0 & 0 & 1 & -1 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C^2}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right). \]

6.5.3.1 Ejemplo 3: Estudio de la influencia de las citoquinas en ratones viejos

El objetivo de este experimento es estudiar el efecto de las citoquinas sobre la estimulación de una sustancia (LPS) y ver como esta relación se ve afectada por la edad.

Se trata de un modelo de un un factor (tratamiento, asignable por el individuo) y un bloque (edad, no asignable, prefijada en cada ratón). En la práctica el modelo equivale al del análisis de la varianza de dos factores.

El diseño experimental para este estudio fue el siguiente:

Array	Tratamiento	Edad
1	LPS	Viejo
2	LPS	Joven
3	Medio	Viejo
4	Medio	Joven
5	LPS	Viejo
6	LPS	Joven
7	Medio	Viejo
8	Medio	Joven
9	LPS	Viejo
10	LPS	Joven
11	Medio	Viejo
12	Medio	Joven

Se utilizáron microarrays de un color (Affymetrix).
Cada condición se replicó tres veces.
Las preguntas específicas a responder:
- ¿Cual es el efecto del tratamiento en ratones viejos?
- ¿Cual es el efecto del tratamiento en ratones jovenes?
- ¿En que genes el efecto es diferente?.

En este caso se pueden considerar distintos modelo lineales derivados del hecho de que este experimento admite diferente parametrizaciones:

Factores separados con 2 niveles cada uno para tratamiento (LPS/Med), Edad (Joven/Viejo) y su interacción: \[Y_{ijk}=\underbrace{\alpha_{i}}_{\text{Trat}}+\underbrace{\beta_{j}}_{\text{Edad}}+\underbrace{\gamma_{ij}}_{\text{Interacción}}+\varepsilon_{ijk},\, i=1,2,\, j=1,2,\, k=1,2,3\] Esta primera parametrización parece natural pero es más complicado confiar en ella para responder las preguntas propuestas.
Un factor combinado con 4 niveles
(LPS.Aged, Med.Aged, LPS.Young, Med.Young)

\[ Y_{ij}=\alpha{i} +\varepsilon_{ij}, \quad i=1,...,4,\, j=1,2,3 . \] Esta parametrización parece más rígida pero se adapta mejor para responder a las preguntas planteadas.

Aquí se adopta la segunda parametrización.

Parámetros del modelo: \[ \begin{array}{ccc} \alpha_1&=& \mathbf{E} (log{LPS.Aged} ),\ \alpha_2= \mathbf{E} (log{Med.Aged} ),\\ \alpha_3&=&\mathbf{E} (log{LPS.Young} ),\, \alpha_4=\mathbf{E} (log{Med.Young} ). \end{array} \]
Contrastes: Preguntas que interesa responder: \[ \begin{array}{ccc} \beta_1^1 &=& \alpha_3-\alpha_1,\quad \mbox{Efecto del tratamiento en ratones viejos} \\ \beta_2^1 &=& \alpha_4-\alpha_2,\quad \mbox{Efecto del tratamiento en ratones jóvenes} \\ \beta_3^1 &=& (\alpha_3-\alpha_1)- (\alpha_2-\alpha_4),\quad \mbox{Interacción: diferencia entre efectos} \end{array} \]
Modelo lineal:

Modelo, $\mathbf{y=X\alpha+\varepsilon}$, y contrastes $\mathbf{C'\beta}$

\[ \left( \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \\ y_7 \\ y_8 \\ y_9 \\ y_{10} \\ y_{11} \\ y_{12} \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \end{array} \right)}_{\text{Matriz de Diseño}, \mathbf{X}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right) + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6\\ \varepsilon_7\\ \varepsilon_8 \end{array} \right) \]

\[ \left( \begin{array}{r} \beta_1^1 \\ \beta_2^1 \\ \beta_3^1 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} -1 & 0 & 1 & 0\\ 0 & -1 & 0 & 1\\ -1 & 1 & 1 & -1 \end{array} \right)}_{\text{Matriz de Contraste}, \mathbf{C}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right) \]

6.5.4 Estimación e inferencia con el modelo lineal

Una vez se ha expresado el experimento como un modelo lineal: \[ \mathbf{E(y_g)=X\alpha_g},\quad \text{var}(y_g)=W_g\sigma_g, \] es posible usar la teoría estándar del modelo lineal (ver (faraway:2004?)) para obtener:

Estimación de los parámetros: $\hat \alpha_g( \approx \mathbf{\alpha})$.
Desviación estándar de las estimaciones: $\hat \sigma_g=s_g( \approx \sigma)$.
Error estándar de las estimaciones:$\widehat {\text{var}\hat \alpha_g}=V_g\,s_g^2$. Estas estimaciones son la base para realizar inferencia sobre $\alpha$ i.e. test $H_0:\ \alpha =0?$, basado en el hecho que: \[ t_{gj}=\frac{\alpha_{gj}}{s_g\sqrt{v_{gj}}}\sim \text{Distribución de Student}. \] De forma análoga pueden derivarse fórmulas para $\alpha_1-\alpha_2$, es decir para decidir acerca de las comparaciones.

Los procedimientos de estimación y de inferencia no dependen de qué parametrización se ha adoptado, a pesar de que distintas parametrizaciones piedan dar lugar a distintos valores numéricos..

6.5.4.1 Fortaleza y debilidades del modelo lineal

El enfoque del modelo lineal es flexible y potente:

Se puede adaptar a situaciones diferentes y complejas.
Siempre produce buenas estimaciones (“BLUE”).
Si las suposiciones son ciertas proporciona una base para la inferencia.

Por otro lado hay que tener en cuenta que, si las suposiciones no se cumplen, entonces las conclusiones deben tomarse con precaución.

Y lo que es peor, aún siendo válidas las suposiciones del modelo los resultados pueden verse afectados por el tamaño de la muestra de forma que, en muestras pequeñas donde pueden haber variaciones más grandes es fácil que se obtengan valores $t$ no significativos o, por el contrario, excesivamente significativos (si la variación es muy reducida).

La metodología desarrollada por Smyth (Smyth, Michaud, and Scott (2005)) basada en los resultados de L”onsted & Speed (Speed (2003)) está dirigida a abordar cómo hacer frente a estas debilidades.

6.5.5 Modelos lineales para Microarrays

Smyth (Smyth, Michaud, and Scott (2005)) considera el problema de la identificación de genes que se expresan diferencialmente en las condiciones especificadas en el diseño de experimentos de microarrays. Como hemos dicho repetídamente la variabilidad de los valores de expresión difiere entre genes, pero la naturaleza paralela de la inferencia en microarrays sugiere la posibilidad de usar la información de todos los genes a la vez para mejorar la estimación de los parámetros, lo que puede llevar a una inferencia más fiable.

Básicamente lo que propone Smyth (Smyth, Michaud, and Scott (2005)) es una solución en tres pasos:

Se plantea el problema como un model lineal con una componente bayesiana ya que se supone que los mismos parámetros a estimar son variables (no constantes) con distribuciones prior que se estimaran a partir de la información de todos los genes.
A continuación se obtienen las estimaciones de los parámetros del modelo. La aproximación utilizada garantiza que estos estimadores tienen un comportamiento robusto incluso para pequeño número de arrays.
Finalmente se calcula un “odd-ratio” que viene a ser la probabilidad de que un gen esté diferencialmente expresado frente a la de que no lo esté y se asocia este valor denominado estadístico $B$ con un estadístico $t$ moderado y su p–valor. \[ B=\log \frac{P[\text{Afectado}|M_{ij}]}{P[\text{No Afectado}|M_{ij}]}, \] gen=i $(i=1...N)$, réplica=j $(j=1,...,n)$.

El hecho de trabajar con logaritmos permite poner el punto de corte en el cero: A mayor valor positivo más probable es que el gen esté diferencialmente expresado. A mayor valor negativo, más probable es que no lo esté

6.5.6 Implementación y ejemplos

El paquete de Bioconductor limma al que hemos hecho referencia en los párrafos anteriores implementa el método de Smyth para la regularización de la varianza lo que lo ha convertido en muy popular entre los usuarios de microarrays

El código siguiente muestra como se crea la matriz del diseño y los contrastes para realizar el análisis mediante modelos lineales;

Empezamos por la matriz de diseño.

##                      Aged.LPS Aged.MED Young.LPS Young.MED
## Aged LPS 80L.CEL            1        0         0         0
## Aged LPS 86L.CEL            1        0         0         0
## Aged LPS 88L.CEL            1        0         0         0
## Aged Medium 81m.CEL         0        1         0         0
## Aged Medium 82m.CEL         0        1         0         0
## Aged Medium 84m.CEL         0        1         0         0
## Young LPS 75L.CEL           0        0         1         0
## Young LPS 76L.CEL           0        0         1         0
## Young LPS 77L.CEL           0        0         1         0
## Young Medium 71m.CEL        0        0         0         1
## Young Medium 72m.CEL        0        0         0         1
## Young Medium 73m.CEL        0        0         0         1
## attr(,"assign")
## [1] 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$lev
## [1] "contr.treatment"

Usando los nombres de las columnas de la matriz de diseño creamos los contrastes

require(limma)
cont.matrix <- makeContrasts (
      LPS.in.AGED=(Aged.LPS-Aged.MED),
      LPS.in.YOUNG=(Young.LPS-Young.MED),
      AGE=(Aged.MED-Young.MED),
      levels=design)
cont.matrix

##            Contrasts
## Levels      LPS.in.AGED LPS.in.YOUNG AGE
##   Aged.LPS            1            0   0
##   Aged.MED           -1            0   1
##   Young.LPS           0            1   0
##   Young.MED           0           -1  -1

Una vez definida la matriz de diseño y los contrastes podemos pasar a estimar el modelo, estimar los contrastes y realizar las pruebas de significación que nos indiquen, para cada gen y cada comparación, si puede considerarse diferencialmente expresado.

La penúltima instrucción ejecuta el proceso de regularización de la varianza utilizando modelos de Bayes empíricos para combinar la información de toda la matriz de datos y de cada gen individual y obtener estimaciones de error mejoradas.

El análisis proporciona los estadísticos de test habituales como Fold–change $t$-moderados o $p$-valores ajustados que se utilizan para ordenar los genes de mas a menos diferencialmente expresados.

A fin de controlar el porcentaje de falsos positivos que puedan resultar del alto numero de contrastes realizados simultaneamente los p–valores se ajustan de forma que tengamos control sobre la tasa de falsos positivos utilizando el metodo de Benjamini y Hochberg (Benjamini and Hochberg (1995)).

La funcion topTable genera para cada contraste una lista de genes ordenados de mas a menos diferencialmente expresados.

topTab_LPS.in.AGED <- topTable (fit.main, number=nrow(fit.main), coef="LPS.in.AGED", adjust="fdr")
topTab_LPS.in.YOUNG <- topTable (fit.main, number=nrow(fit.main), coef="LPS.in.YOUNG", adjust="fdr")
topTab_AGE  <- topTable (fit.main, number=nrow(fit.main) , coef="AGE", adjust="fdr")

Se puede visualizar los resultados mediante un volcano plot que, en este caso, representa en abscisas los cambios de expresión en escala logarítmica y en ordenadas el estadístico $B$ en vez de el “menos logaritmo” del p-valor.

coefnum = 1
opt <- par(cex.lab = 0.7)
volcanoplot(fit.main, coef=coefnum, highlight=10, names=fit.main$ID,
            main=paste("Differentially expressed genes",colnames(cont.matrix)[coefnum], sep="\n"))
abline(v=c(-1,1))

par(opt)

7 Caso “Estrógeno”: Selección de genes diferencialmente expresados asociados con la resistencia a estrógenos.

7.1 Introducción

El análisis de datos de microarrays suele proceder, en general secuencialmente, a través de una serie de etapas tal y como se comentado en capítulos anteriores.

En cada una de las etapas pueden llevarse a cabo diversos procesos, con diversas variantes de métodos parecidos.A lo largo de la primera década del 2000 se ha desarrollado el proyecto Bioconductor que ha impulsado la creación de cientos de paquetes de R que implementan casi todas las capacidades posibles de análisis de microarrays. Aunque nadie puede conocer todos los paquetes con tan sólo conocer algunos de ellos es posible llevar a cabo potentes análisis con relativa facilidad.

El objetivo de este capítulo es ilustrar de forma breve pero completa como se hace un análisis de microarrays usando las facilidades que ofrecen Bioconductor y R.

Como se ha dicho, uno de los problemas en este tipo de estudios es que en cada etapa se puede proceder de varias formas lo que da lugar a un asfixiante número de posibilidades especialmente para el neófito. Con el fin de evitar este problema de momento se describe una sola de las opciones posibles para cada paso, lo que constituye un proceso “al estilo Bioconductor”

7.2 El estudio “Estrógeno”

El ejemplo analizado consiste en unos datos disponibles en forma de paquete (estrogen) en Bioconductor. Se trata de un estudio en el que se estudia la influencia del tratamiento con estrógeno y del tiempo transcurrido desde el tratamiento en la expresión de genes asociados con cancer de mama. En capítulos anteriores se han utilizado fragmentos de este ejemplo descrito en el capítulo here en el epígrafe corresponiente a los ejemplos (ejemplo here).

El diseño experimental para este estudio consiste en un diseño factorial de 2 factores “Estrogeno” (Presente o Ausente), “Tiempo” (10h o 48h). La documentación del paquete estrogen ofrece más detalles acerca del diseño y las variables utilizadas.

Si el paquete no se encuentra instalado puede instalarse con la instrucción BiocManager::install que se descarga de la web de Bioconductor.

if (!(require("estrogen", character.only=T))){
    BiocManager::install("estrogen")
    }

En lo que sigue supondremos que se ha llevado a cabo una instalación estándar de Bioconductor es decir se han ejecutado las instrucciones:

BiocManager::install()

7.2.1 Directorios y opciones de trabajo

Para facilitar supondremos que trabajamos en un directorio escogido por nosotros y cuya localización se asigna a la variable workingDir. Los datos se copiaran en un subdirectorio del anterior denominado “data” que se almacenará en la variable dataDir y los resultados se almacenarán en un directorio “results” cuyo nombre completo se almacenará en la variable resultsDir.

workingDir <-getwd()
if (!file.exists("datos")) system("mkdir datos")
if (!file.exists("datos/estrogen")) system("mkdir datos/estrogen")
if (!file.exists("results")) system("mkdir results")
dataDir <-file.path(workingDir, "datos/estrogen")
resultsDir <- file.path(workingDir, "results")
setwd(workingDir)

options(width=80)
options(digits=5)

7.3 Obtención y lectura de los datos

Para un análisis de datos de microarrays de Affymetrix se necesitan los archivos de imágenes escaneadas (“.CEL”) y un archivo en el que se asigne una condición experimental a cada archivo.

7.3.1 Los datos

Una ventaja de este ejemplo es que los datos se encuentran disponibles tras instalar el paquete estrogen por lo que se pueden copiar un directorio de trabajo o simplemente trabajar con ellos en el directorio original. El directorio en que se encuentran los archivos .CEL es:

library(estrogen)
estrogenDir <- system.file("extdata", package = "estrogen")
print(estrogenDir)

## [1] "/home/alex/R/x86_64-pc-linux-gnu-library/4.1/estrogen/extdata"

Para realizar el análisis es preciso copiar todos los archivos del directorio “extdata” a nuestro directorio de trabajo “data.”

ATENCION: Esta operación que depende de la instalación y del sistema operativo. Por ejemplo para copiarlos desde R en linux se escribirá:

system(paste ("cp ", estrogenDir,"/* ./data", sep=""))

Los archivos copiados son

Los ocho archivos .CEL para el análisis
Un archivo defectuoso para que se vea como debe salir un “array malo” en el control de calidad.
Un archivo de covariables o “targets” que se usará al leer los archivos .CEL incorporar la información de las covariables en el objeto que se creará.

7.3.2 Lectura de los datos

El proceso de leer los datos puede parecer un poco extraño a primera vista pero la idea es simple:

En primer lugar creamos algunas estructuras de datos que contienen la información sobre las variables y - opcionalmente - sobre las anotaciones y el experimento y

A continuación nos basamos en estas estructuras para leer los datos y crear los objetos principales que se utilizarán para el análisis.

library(Biobase)
library(affy)
sampleInfo <- read.AnnotatedDataFrame(file.path(estrogenDir,"targLimma.txt"),
    header = TRUE, row.names = 1, sep="\t")
fileNames <- pData(sampleInfo)$FileName
rawData <- read.affybatch(filenames=file.path(estrogenDir,fileNames),
                          phenoData=sampleInfo)

En este ejemplo se incluye un archivo defectuoso con fines didácticos, llamado “badcel.” Para ilustrar las diferencias entre archivos correctos e incorrectos se puede leer en otro objeto.

require(affy)
sampleInfoWithBad <- read.AnnotatedDataFrame(file.path(dataDir,"phenoDataWithBad.txt"),
    header = TRUE, row.names = NULL, sep="\t")
fileNames <- pData(sampleInfoWithBad)$FileName
rawData.wrong <- read.affybatch(filenames=file.path(dataDir,fileNames),
                          phenoData=sampleInfoWithBad)

7.4 Exploración, Control de Calidad y Normalización

Tras leer los datos pasamos al preprocesado. Aunque puede interpretarse de distintas formas esta fase suele consistir en

Realizar algunos gráficos con los datos para hacerse una idea de como ha resultado el experimento.
Realizar un control de calidad más formal.
Normalizar y, en el caso de Affymetrix, resumir las expresiones

7.4.1 Exploración y visualización

La exploración previa puede hacerse paso a paso o en bloque si se utiliza algún paquete como arratQualityMetrics. En la ayuda de este paquete se encuentra una descripción de los análisis básicos que permite realizar.

El primer gráfico que suele hacerse és algun tipo de “densidad” que muestre la distribución de las señales en los datos “crudos,” es decir sin normalizar. Estos gráficos también permiten hacerse una idea de si las distribuciones de los distintos arrays son similares en forma y posición.

El gráfico de degradación –que no aparece en este caso, ya que daba problemeas al representarlo– permite hacerse una idea de como ha sido el proceso de hibridación de las muestras. Lineas paralelas sugieren una calidad similar.

##        low10A low10B  hi10A   hi10B   low48A   low48B   hi48A    hi48B
## slope  -0.103 -0.217 -0.129 -0.3810 -0.50400 -0.55300 -0.3510 -0.66700
## pvalue  0.550  0.194  0.394  0.0341  0.00482  0.00109  0.0464  0.00032

EL código para el gráfico de degradación sería:

El diagrama de cajas muestra, como el histograma, una idea de la distribución de los datos.

Finalmente un cluster jerárquico seguido de un dendrograma nos puede ayudar a hacernos una idea de si las muestras se agrupan por condiciones experimentales.

Si lo hacen es bueno, pero si no, no es necesariamente indicador de problemas, puesto que es un gráfico basado en todo los datos.

7.4.2 Control de calidad

Las exploraciones anteriores nos proporcionan una idea de como son los datos.

Se pueden realizar controles de calidad más estrictos como:

Los controles de calidad estándar de Affymetrix, descritos en el paquete arrayQualityMetrics.

El paquete arrayQualityMetrics encapsula unos cuantos análisis, de forma que con una instrucción se pueden realizar todos los análisis y enviar la salida a un directorio.

if(!require(arrayQualityMetrics)) BiocManager::install("arrayQualityMetrics")
library(arrayQualityMetrics)
arrayQualityMetrics(rawData, outdir = "Informe_de_calidad_para_los_datos_del_caso_ESTROGENO")
arrayQualityMetrics(rawDataBad, outdir = "Informe_de_calidad_(2)_para_los_datos_del_caso_ESTROGENO")

Estos informes contienen gráficos, y explicaciones de como interpretarlos. Los gráficos llevan también indicadores para determinar si una muestra determinada puede tener algún tipo de problema.

En el código de ejemplo hemos realizado dos veces el análisis, una con los datos “buenos” y otra incluyendo el array defectuoso (“bad.cel”). Como puede verse en las tablas resumen, y en los directorios de resultados el array defectuoso es detectado por la mayoría de pruebas.

Para concluir con esta sección merece la pena recordar que todos estos gráficos son exploratorios. Raramente se descarta un array si tan sólo un grafico lo sugiere.

7.4.3 Normalizacion y Filtraje

Una vez realizado el control de calidad se procede a normalizar los datos y sumarizarlos.

Hecho esto puede realizarse un filtraje no específico con el fin de eliminar genes que constituyen básicamente “ruído,” bien porque sus señales son muy bajas o bien porque apenas varían entre condiciones, por lo que no aportan nada a la selección de genes diferencialmente expresados.

La normalización puiede hacerse por distintos métodos (MAS5, VSN, RMA, GCRMA, …) En este ejemplo se utilizará el método RMA pero no se realizará filtraje alguno. Esto puede implicar quizás que para seleccionar genes diferencialmente expresados basándose en el ajuste de p-valores debamos utilizar criterios menos restrictivos que si hubiéramos filtrado, pero tiene la ventaja de eliminar un paso que, en el mejor de los casos, resulta controvertido.

El procesado mediante RMA implica un proceso en tres etapas:

Corrección de fondo (el RMA hace precisamente esto).
Normalización para hacer los valores de los arrays comparables.
Summarización de las diversas sondas asociadas a cada grupo de sondas para dar un único valor.

require(affy)
eset_rma <- rma(rawData)

## Background correcting
## Normalizing
## Calculating Expression

Si se desea normalizar únicamente si los datos normalizados no estan disponibles puede utilizarse el código siguiente:

require(affy)
normalize <- T
if(normalize){
  eset_rma <- rma(rawData)
  save(eset_rma, file=file.path(dataDir,"estrogen-normalized.Rda"))
}else{
  load (file=file.path(dataDir,"estrogen-normalized.Rda"))
}

La normalización hace que los valores de los arrays sean comparables entre ellos, aunque los distintos métodos situan los valores en escalas distintas, por lo que lo que no resulta directamente comparable son los valores normalizados por distintos métodos.

Un boxplot de los valores normalizados muestra que los valores ya están en una escala en donde se pueden comparar.

Ahora bien, esto no debe considerarse sinónimo de buena calidad porque si se aplica la normalización por cuantiles el gráfico aparecerá así independientemente de otros factores que puedan afectar la calidad.

boxplot(eset_rma,main="Boxplot for RMA-normalized expression values ",
        names=sampleNames, cex.axis=0.7, col=info$grupo+1,las=2)

El código siguiente, que no se ejecuta muestra como se podrían aplicar distintos métodos para lego compararlos entre ellos.

eset_mas5 <- mas5(rawData)  # Uses expresso (MAS 5.0 method) much slower than RMA!
stopifnot(require(gcrma))
eset_gcrma <- gcrma(rawData) # The 'library(gcrma)' needs to be loaded first.
stopifnot(require(plier))
eset_plier <- justPlier(rawData, normalize=T) # The 'library(plier)' needs to be loaded first.
compara <-data.frame(RMA=exprs(eset_rma)[,1], MAS5 =exprs(eset_mas5)[,1],
                    GCRMA=exprs(eset_gcrma)[,1], PLIER =exprs(eset_plier)[,1])
pairs(compara)

7.4.3.1 Filtraje

El filtraje no específico permite eliminar los genes que varían poco entre condiciones o que deseamos quitar por otras razones como por ejemplo que no disponemos de anotación para ellos. La función nsFilter permite eliminar los genes que, o bien varían poco, o bien no se dispone de anotación para ellos.

Si al filtrar deseamos usar las anotaciones, o la falta de ellas, como criterio de filtraje debemos disponer del correspondiente paquete de anotaciones.

require(genefilter)
filtered <- nsFilter(eset_rma, require.entrez=TRUE,
         remove.dupEntrez=TRUE, var.func=IQR,
         var.cutoff=0.5, var.filter=TRUE,
         filterByQuantile=TRUE, feature.exclude="^AFFX")

La función nsFilter devuelve los valores filtrados en un objeto expressionSet y un informe de los resultados del filtraje.

names(filtered)

## [1] "eset"       "filter.log"

class(filtered$eset)

## [1] "ExpressionSet"
## attr(,"package")
## [1] "Biobase"

print(filtered$filter.log)

## $numDupsRemoved
## [1] 2888
## 
## $numLowVar
## [1] 4388
## 
## $numRemoved.ENTREZID
## [1] 942
## 
## $feature.exclude
## [1] 19

eset_filtered <-filtered$eset

Podemos grabar el objeto eset_rma y los datos filtrados para su posterior uso.

save(eset_rma, eset_filtered, file=file.path(resultsDir, "estrogen-normalized.Rda"))

Despues del filtraje han quedado 4409 genes disponibles para analizar.

7.5 Selección de genes diferencialmente expresados

Como en las etapas anteriores la selección de genes diferencialmente expresados (GDE) puede basarse en distintas aproximaciones, desde la $t$ de Student al programa SAM pasando por multitud de variantes.

En este ejemplo se aplicará la aproximación presentada por Smyth et al. (2004) basado en la utilización del modelo lineal general combinada con un método para obtener una estimación mejorada de la varianza.

7.5.1 Análisis basado en modelos lineales

El capítulo 6 de tse manual o el manual del programa limma contienen explicaciones detalladas sobre construir un modelo lineal para este problema y como utilizarlo para seleccionar genes diferencialmente expresados.

7.5.1.1 Matriz de diseño

El primer paso para el análisis es crear la matriz de diseño.

La situación discutida en este ejemplo se puede modelizar de dos formas, tal como se discute en la presentación citada más arriba:

Como un modelo de dos factores Estrogeno(Pres/Aus) Tiempo (10h/48h) con interacción.
Como un modelo de un factor con cuatro niveles (Pres.10h/Pres.48h/Aus.10h/Aus.48h).

Tal como se describe en el manual de limma la segunda parametrización resulta a menudo más comoda, a pesar de parecer menos intuitiva, porque permite formular con más facilidad que la de dos factores las preguntas que típicamente interesan a los investigadores.

El modelo lineal para este estudio será:

\[ \left( \begin{array}{r} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \\ y_7 \\ y_8 \end{array} \right) = \underbrace{\left( \begin{array}{rrrr} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \end{array} \right)}_{\mbox{Design Matrix}, \mathbf{X}} \left( \begin{array}{c} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{array} \right) + \left( \begin{array}{r} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3\\ \varepsilon_4 \\ \varepsilon_5\\ \varepsilon_6\\ \varepsilon_7\\ \varepsilon_8 \end{array} \right) \]

Los parámetros del modelos representan las cuatro combinaciones tiempo/estrogeno.

\[\begin{array}{ccc} \alpha_1&=& \mathbf{E} (log{Abs.10h} ),\\ \alpha_2&=& \mathbf{E} (log{Pres.10h} ),\\ \alpha_3&=& \mathbf{E} (log{Abs.48h} ),\\ \alpha_4&=& \mathbf{E} (log{Pres.48h} ). \end{array}\]

La matriz de diseño puede definirse manualmente o a partir de un factor creado específicamente para ello.

Manualmente, seria:

design.1<-matrix(
c(1,1,0,0,0,0,0,0,
  0,0,1,1,0,0,0,0,
  0,0,0,0,1,1,0,0,
  0,0,0,0,0,0,1,1),
nrow=8,
byrow=F)
colnames(design.1)<-c("neg10h", "est10h", "neg48h", "est48h")
rownames(design.1) <-  c("low10A", "low10B", "hi10A" , "hi10B",  "low48A", "low48B", "hi48A" , "hi48B")
print(design.1)

##        neg10h est10h neg48h est48h
## low10A      1      0      0      0
## low10B      1      0      0      0
## hi10A       0      1      0      0
## hi10B       0      1      0      0
## low48A      0      0      1      0
## low48B      0      0      1      0
## hi48A       0      0      0      1
## hi48B       0      0      0      1

Alternativamente puede crearse la matriz de diseño a partir de la información sobre las condiciones contenida en el phenoData, siempre que exista un campo adecuado para ello.

En este caso la columna Target se ha creado para utilizarla con esta finalidad. El objeto phenoData puede recrearse a partir del archivo original o extrayéndolo del objeto ExpresionSet que contiene los datos y las covariables.

##        neg10h est10h neg48h est48h
## low10A      1      0      0      0
## low10B      1      0      0      0
## hi10A       0      1      0      0
## hi10B       0      1      0      0
## low48A      0      0      1      0
## low48B      0      0      1      0
## hi48A       0      0      0      1
## hi48B       0      0      0      1
## attr(,"assign")
## [1] 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$lev
## [1] "contr.treatment"

Ambas matrices, design, design.1 resultan iguales.

7.5.1.2 Contrastes

Dado un modelo lineal definido a traves de una matriz de diseño pueden formularse las preguntas de interes como contrastes es decir comparaciones entre los parametros del modelo.

Cada parametrizacion distinta requerira de unos contrastes diferentes para las mismas preguntas, por lo que habitualmente se utilizara la parametrizacion que permita formular de forma mas clara las comparaciones de interes.

En este caso interesa estudiar

Efecto del estrogeno al inicio del tratamiento
Efecto del estrogeno al cabo de un tiempo del tratamiento
Efecto del tiempo en ausencia de estrogeno

Esto se puede formular facilmente con la parametrizacion adoptada. \[\begin{array}{ccc} \beta_1^1 &=& \alpha_2-\alpha_1,\quad \mbox{Efecto del estrogeno pasadas 10 horas} \\ \beta_2^1 &=& \alpha_4-\alpha_3,\quad \mbox{Efecto del estrogeno pasadas 48 horas} \\ \beta_3^1 &=& \alpha_3-\alpha_1,\quad \mbox{Efecto del tiempo en ausencia de Estrogeno} \end{array}\]

cont.matrix <- makeContrasts (
      Estro10=(est10h-neg10h),
      Estro48=(est48h-neg48h),
      Tiempo=(neg48h-neg10h),
      levels=design)
cont.matrix

##         Contrasts
## Levels   Estro10 Estro48 Tiempo
##   neg10h      -1       0     -1
##   est10h       1       0      0
##   neg48h       0      -1      1
##   est48h       0       1      0

7.5.1.3 Estimación del modelo y selección de genes

require(limma)
fit<-lmFit(eset_rma, design)
fit.main<-contrasts.fit(fit, cont.matrix)
fit.main<-eBayes(fit.main)

El método implementado en amplía el análisis tradicional utilizando modelos de Bayes empíricos para combinar la información de toda la matriz de datos y de cada gen individual y obtener estimaciones de error mejoradas.

La funcion topTable genera para cada contraste una lista de genes ordenados de mas a menos diferencialmente expresados.

topTabEstro10 <- topTable (fit.main, number=nrow(fit.main), coef="Estro10", adjust="fdr")
topTabEstro48 <- topTable (fit.main, number=nrow(fit.main), coef="Estro48", adjust="fdr")
topTabTiempo  <- topTable (fit.main, number=nrow(fit.main) , coef="Tiempo", adjust="fdr")

Una forma de visualizar los resultados es mediante un volcano plot que representa en abscisas los cambios de expresión en escala logarítmica y en ordenadas el “menos logaritmo” del p-valor o alternativamente el estadístico $B$.

7.5.2 Comparaciones múltiples

Cuando se realizan varias comparaciones a la vez puede resultar importante ver que genes cambian simultáneamente en más de una comparación. Si el número de comparaciones es alto también puede ser necesario realizar un ajuste de p-valores entre las comparaciones, distinto del realizado entre genes.

La función decidetests permite realizar ambas cosas. En este ejemplo no se ajustaran los p–valores entre comparaciones. Tan solo se seleccionaran los genes que cambian en una o más condiciones.

EL resultado del análisis es una tabla res que para cada gen y cada comparación contiene un 1 (si el gen esta sobre-expresado o “up” en esta condicion), un 0 (si no hay cambio significativo) o un -1 (si esta “down”-regulado).

Para resumir dicho análisis podemos contar qué filas tienen como mínimo una celda distinta de cero:

sum.res.rows<-apply(abs(res),1,sum)
res.selected<-res[sum.res.rows!=0,]
print(summary(res))

##        Estro10 Estro48 Tiempo
## Down        23      99     27
## NotSig   12512   12375  12591
## Up          90     151      7

Un diagrama de Venn permite visualizar la tabla anterior sin diferenciar entre genes “up” o “down” regulados.

\begin{figure}

vennDiagram (res.selected[,1:3], main="Genes in common #1", cex=0.9)

\end{figure}

8 Después de la selección: Análisis de listas de genes

8.1 Introducción

El resultado de un análisis de microarrays como los descritos hasta el momento suele ser una “lista de genes,” es decir una archivo o una tabla con los identificadores de los genes considerados diferencialmente expresados en una o mÃ¡s comparaciones.

La pregunta obvia que los investigadores se plantean al disponer de estas listas es ¿cuál es su significado biológico?.

Para el bioinfomatico eta pregunta se replantea de la siguiente manera: ¿Que se puede hacer para convertir -de forma más o menos automática– la lista obtenida en algún tipo de información biológicamente relevante. Esta cuestión ha llevado al desarrollo de métodos y herramientas que de forma general se agrupan bajo el título de métodos para el análisis de listas de genes o métodos de análisis de la significación biológica.

Siendo este un campo muy extenso no lo vamos a desarrollar en profundidad sinó que nos limitaremos a considerar brevemente algunas aproximaciones sencillas que pueden contribuir a una mejor interpretación de los resultados de un experimento de microarrays.

Los aspectos que trataremos serán:

Análisis comparativo de listas de genes.
Anotación de los resultados.
Visualización de la matriz de expresión para los genes seleccionados en una o más comparaciones.
Análisis básico de la significación biológica: Gene Enrichment Analysis

8.2 Análisis comparativo de listas de genes

Muchos estudios analizan el comportamiento de los genes bajo varios tratamientos o condiciones experimentales.

En estos casos puede ser interesante ver como cambia un gen bajo distintos tratamientos o bien ver qué genes son afectados de forma parecida o distinta bajo los mismos. La forma habitual de hacer esta comparación es mediante algún programa que permita comparar los elementos de dos o más listas. En Bioconductor, por ejemplo el paquete limma tiene algunas funciones que permiten seleccionar los genes cambiados simultaneamente en dos o más condiciones.

8.2.1 Ejemplo: comparación de listas en el caso `celltypes`

Supondremos que tras las comparaciones del capítulo anterior disponemos de los valores de expresión y las listas de genes obtenidas de las tres comparaciones del caso cellypes. Ambos están almacenados en archivos binarios que podemos recuperar con la instrucción load.

require(Biobase)
require(limma)
load("datos/celltypes/celltypes-normalized.rma.Rda")
load("datos/celltypes/celltypes-fit.main.Rda")
topTab_LPS.in.AGED <- topTable (fit.main, number=nrow(fit.main), coef="LPS.in.AGED", adjust="fdr",lfc=2)
topTab_LPS.in.YOUNG <- topTable (fit.main, number=nrow(fit.main), coef="LPS.in.YOUNG", adjust="fdr",lfc=2)
topTab_AGE  <- topTable (fit.main, number=nrow(fit.main) , coef="AGE", adjust="fdr", lfc=2)

La función decidetests permite realizar ambas cosas. En este ejemplo no se ajustaran los p–valores entre comparaciones. Tan solo se seleccionaran los genes que cambian en una o más condiciones.

Para resumir dicho análisis podemos contar qué filas tienen como mínimo una celda distinta de cero:

sum.res.rows<-apply(abs(res),1,sum)
res.selected<-res[sum.res.rows!=0,]
print(summary(res))

##        LPS.in.AGED LPS.in.YOUNG  AGE
## Down           556          561  558
## NotSig        9281         9239 8740
## Up             417          454  956

En vista de estos valores podemos aplicar otros criterios de selección, por ejemplo genes con un p-valor ajustado inferior a 0.0001 y log Fold change mayor o igual a 2._ Este criterio combina a la vez la significación estadística y la significación biológica por lo que, en un estudio real sería probablemente el escogido_.

##        LPS.in.AGED LPS.in.YOUNG   AGE
## Down            71           83    45
## NotSig       10106        10083  9979
## Up              77           88   230

Un diagrama de Venn permite visualizar la tabla anterior sin diferenciar entre genes “up” o “down” regulados.

vennDiagram (res.selected[,1:3], main="Genes in common #1", cex=0.9)

Figure 8.1: Número de genes seleccionado en cada comparación

8.3 Anotación de los resultados

La identificación de los genes seleccionados puede resultar más sencilla para el especialista en un campo si se utilizan nombres estándar como el simbolo del gen o “gene symbol.” Ahora bien hemos de recordar que los microarrays suele ser un producto comercial en los que la empresa que los crea decide que sondas de cada gen o exón coloca en cada posición del array. Esto significa que para poder indicar el nombre de cada gen seleccionado es preciso disponer de algún tipo de tabla de anotaciones que relacione cada sonda con los genes a los que apunta. Esta tabla debería ser proporcionada por las companñias que producen los arrays pero es algo que queda a su criterio y, en el mejor de los casos, cada empresa puede seguir unos criterios distintos.

Afortunadamente si se trabaja con Bioconductor existe una solución homogenea a este problema. Se trata de los paquetes de anotaciones creados y mantenidos por el equipo de Bioconductor a partir de la información publicada por cada compañia pero con una forma común de acceso. Es decir aunque las distintas empresas utilicen fotmatos distintos desde el punto de vista del usuario de Bioconductor esto no importará: si existe el paquete se utiliza como todos los demás.

De hecho Bioconductor incorpora varios tipos de paquetes de anotaciones:

Paquetes de anotaciones de microarrays: Asocian a vada sonda la información disponible para ella en distintas bases de datos. Consisten en una serie de tablas SQL (de hecho cada paquete es una base de datos SQL) que asocian las sondas de cada gen con sus correspondencias en distintas bases de datos.
Paquetes de anotaciones de organismos: Contienen asociaciones similares a las anteriores pero las claves principales de cada tabla no son las sondas de un microarray determinado sino los genes de un organismo dado. Es decir en estos paquetes cada tabla asocia el identificador entrez del gen con sus correspondientes identificadores en otras bases de datos.
Paquetes de anotaciones de bases de datos Estos paquetes difieren de los anteriores porque lo que hacen es almacenar en el mismo formato que los anteriores –bases de datos SQL consultables en R– una copia de ciertas bases de datos como la Gene Ontology (GO) o la Kyoto Encyclopedia of Genes and Genomes (KEGG).

Para saber que anotaciones estan disponibles debe cargarse el paquete y llamar la función del mismo nombre. Por ejempl para los microarrays del caso celltypes:

Para saber que anotaciones estan disponibles debe cargarse el paquete y llamar la función del mismo nombre sin el sufijo “.db.”

library(mouse4302.db)
anotData <- capture.output(mouse4302())
print(anotData)

##  [1] "Quality control information for mouse4302:"                  
##  [2] ""                                                            
##  [3] ""                                                            
##  [4] "This package has the following mappings:"                    
##  [5] ""                                                            
##  [6] "mouse4302ACCNUM has 45101 mapped keys (of 45101 keys)"       
##  [7] "mouse4302ALIAS2PROBE has 81336 mapped keys (of 156785 keys)" 
##  [8] "mouse4302CHR has 37432 mapped keys (of 45101 keys)"          
##  [9] "mouse4302CHRLENGTHS has 66 mapped keys (of 66 keys)"         
## [10] "mouse4302CHRLOC has 34742 mapped keys (of 45101 keys)"       
## [11] "mouse4302CHRLOCEND has 34742 mapped keys (of 45101 keys)"    
## [12] "mouse4302ENSEMBL has 34432 mapped keys (of 45101 keys)"      
## [13] "mouse4302ENSEMBL2PROBE has 18067 mapped keys (of 33358 keys)"
## [14] "mouse4302ENTREZID has 37466 mapped keys (of 45101 keys)"     
## [15] "mouse4302ENZYME has 4480 mapped keys (of 45101 keys)"        
## [16] "mouse4302ENZYME2PROBE has 954 mapped keys (of 971 keys)"     
## [17] "mouse4302GENENAME has 37466 mapped keys (of 45101 keys)"     
## [18] "mouse4302GO has 34021 mapped keys (of 45101 keys)"           
## [19] "mouse4302GO2ALLPROBES has 22446 mapped keys (of 22649 keys)" 
## [20] "mouse4302GO2PROBE has 18283 mapped keys (of 18511 keys)"     
## [21] "mouse4302MGI has 37455 mapped keys (of 45101 keys)"          
## [22] "mouse4302MGI2PROBE has 20485 mapped keys (of 70076 keys)"    
## [23] "mouse4302PATH has 10731 mapped keys (of 45101 keys)"         
## [24] "mouse4302PATH2PROBE has 225 mapped keys (of 225 keys)"       
## [25] "mouse4302PMID has 37302 mapped keys (of 45101 keys)"         
## [26] "mouse4302PMID2PROBE has 327939 mapped keys (of 340004 keys)" 
## [27] "mouse4302REFSEQ has 35206 mapped keys (of 45101 keys)"       
## [28] "mouse4302SYMBOL has 37466 mapped keys (of 45101 keys)"       
## [29] "mouse4302UNIPROT has 33001 mapped keys (of 45101 keys)"      
## [30] ""                                                            
## [31] ""                                                            
## [32] "Additional Information about this package:"                  
## [33] ""                                                            
## [34] "DB schema: MOUSECHIP_DB"                                     
## [35] "DB schema version: 2.1"                                      
## [36] "Organism: Mus musculus"                                      
## [37] "Date for NCBI data: 2021-Apr14"                              
## [38] "Date for GO data: 2021-02-01"                                
## [39] "Date for KEGG data: 2011-Mar15"                              
## [40] "Date for Golden Path data: 2021-Feb16"                       
## [41] "Date for Ensembl data: 2021-Feb16"

cat ("... output continues until ", length(anotData), " lines.\n")

## ... output continues until  41  lines.

Si en vez del paquetes del microarray usáramos el paquete de organismo, en este caso del ratón.

if (!(require(org.Mm.eg.db))){
        biocLite("org.Mm.eg.db")
      }
require(org.Mm.eg.db)
anotData <- capture.output(org.Mm.eg())
print(anotData[1:15])

##  [1] "Quality control information for org.Mm.eg:"                
##  [2] ""                                                          
##  [3] ""                                                          
##  [4] "This package has the following mappings:"                  
##  [5] ""                                                          
##  [6] "org.Mm.egACCNUM has 46208 mapped keys (of 73083 keys)"     
##  [7] "org.Mm.egACCNUM2EG has 606678 mapped keys (of 606678 keys)"
##  [8] "org.Mm.egALIAS2EG has 156785 mapped keys (of 156785 keys)" 
##  [9] "org.Mm.egCHR has 72225 mapped keys (of 73083 keys)"        
## [10] "org.Mm.egCHRLENGTHS has 66 mapped keys (of 66 keys)"       
## [11] "org.Mm.egCHRLOC has 26193 mapped keys (of 73083 keys)"     
## [12] "org.Mm.egCHRLOCEND has 26193 mapped keys (of 73083 keys)"  
## [13] "org.Mm.egENSEMBL has 33468 mapped keys (of 73083 keys)"    
## [14] "org.Mm.egENSEMBL2EG has 33358 mapped keys (of 33358 keys)" 
## [15] "org.Mm.egENSEMBLPROT has 9437 mapped keys (of 73083 keys)"

cat ("... output continues until ", length(anotData), " lines.\n")

## ... output continues until  48  lines.

Cada tabla de asociación puede consultarse de diversas formas,

Con las funciones get o mget.
Convirtiéndola en una tabla y extrayendo valores
En algunos casos utilizando funciones específicas como getSYMBOL o getEG (por “Entrez Gene”) cuando exitan.

Por ejemplo si tomamos los cinco primeros genes seleccionados en la comparación “LPS.in.AGED”

top5 <-rownames(topTab_LPS.in.AGED)[1:5]
cat("Usando mget\n")

## Usando mget

geneSymbol5.1 <- unlist(mget(top5, mouse4302SYMBOL))
geneSymbol5.1

## 1450826_a_at   1449450_at   1419607_at 1419681_a_at   1419209_at 
##       "Saa3"      "Ptges"        "Tnf"      "Prok2"      "Cxcl1"

cat("Usando toTable\n")

## Usando toTable

genesTable<- toTable(mouse4302SYMBOL)
rownames(genesTable) <-  genesTable$probe_id
genesTable[top5, 2]

## [1] "Saa3"  "Ptges" "Tnf"   "Prok2" "Cxcl1"

cat("Usando getSYMBOL\n")

## Usando getSYMBOL

require(annotate)
geneSymbol5.3 <- getSYMBOL(top5, "mouse4302.db")
geneSymbol5.3

## 1450826_a_at   1449450_at   1419607_at 1419681_a_at   1419209_at 
##       "Saa3"      "Ptges"        "Tnf"      "Prok2"      "Cxcl1"

Bioconductor dispone de algunos paquetes que permiten aprovechar esta funcionalidad anterior para obtener las anotaciones de cada gen y generar una tabla HTML con enlaces a algunas bases de datos.

De forma sencilla es posible obtener tablas con las anotaciones correspondientes a los genes seleccionados. Si se desea ser más ambicioso es posible generar tablas en las que se combinen hiperenlaces a las anotaciones con los resultados de la selección de genes.

El paquete annafy permite de forma muy simple generar una tabla de anotaciones con hiperenlaces a las bases de datos para cada anotación seleccionada.

La instrucción siguiente crearia una tabla con las anotaciones disponibles para los genes seleccionados en la sección de comparaciones múltiples.

require(annaffy)
genesSelected <- rownames(res.selected)
at <- aafTableAnn(genesSelected, "mouse4302.db")
saveHTML (at, file="results/anotations.html",
          "Annotations for selected genes")

8.4 Visualización de los perfiles de expresión

Tras seleccionar los genes diferencialmente expresados podemos visualizar las expresiones de cada gen agrupándolas para destacar los genes que se encuentran up o down regulados simultáneamente constituyendo perfiles de expresión.

Hay distintas formas de visualización pero aquí tan sólo se presenta el uso de mapas de color o Heatmaps cuyos fundamentos se explican en el capítulo dedicado al descubrimiento de clases.

En primer lugar seleccionamos los genes a visualizar: Se toman todos aquellos que han resultado diferencialmente expresados en alguna de las tres comparaciones.

probeNames<-rownames(res)
probeNames.selected<-probeNames[sum.res.rows!=0]
exprs2cluster <-exprs(eset_rma_filtered)[probeNames.selected,]
colnames(exprs2cluster)<- c("OldLPS80L", "OldLPS86L", "OldLPS88L",
                            "OldMED81m", "OldMED82m", "OldMED84m",
                            "YouLPS75L",  "YouLPS76L",  "YouLPS77L",
                            "YouMED71m", "YouMED72m", "YouMED73m")

Para representar el Heatmap tan sólo necesitamos la matriz de datos resultante.

color.map <- function(grupo) {
  switch(grupo,
         "yellow",
         "red",
         "blue",
         "pink")
}
grupColors <- unlist(lapply(pData(eset_rma_filtered)$grupo, color.map))
heatmap(exprs2cluster,
        cexCol=0.8,
        main="Heatmap para las tres comparaciones de 'celltypes'", cex.main=0.8)

Si se desea realizar mapas de color más sofisticados puede utilizarse el paquete Rpackage{gplots} que implementa una version mejorada en la función heatmap.2

require("gplots")
heatmap.2(exprs2cluster,
          col=bluered(75), scale="row",
          ColSideColors=grupColors, key=TRUE, symkey=FALSE,
          density.info="none", trace="none", cexCol=0.8,
          main="Heatmap de las muestras de 'celltypes'",   cex.main=0.6)

8.5 Análisis de significación biológica

En las secciones anteriores se ha visto como encontrar los identificadores de los genes en distintas bases de datos, lo que permite por ejemplo conocer sus nombres comunes (“gene symbol”).

Otra aproximación razonable es estudiar las funciones de los genes buscando sus anotaciones en bases de datos de anotación funcional como la Gene Ontology (GO), o la Kyoto Encyclopedia of Genes and Genomes.

Por ejemplo las anotaciones en GO para los cinco primeros genes de la lista analizada en <a href=``{r anchorLoc('anotExample')“>here serian:

require(annotate)
(top1 <-rownames(topTab_LPS.in.AGED)[1])

## [1] "1450826_a_at"

(geneSymbol1 <- getSYMBOL(top1, "mouse4302.db"))

## 1450826_a_at 
##       "Saa3"

GOAnots1 <- mget(top1, mouse4302GO)
for (i in 1:length(GOAnots1)){
  for (j in 1:length(GOAnots1[[i]])){
    GOAnot <- GOAnots1[[i]][[j]][[1]]
    cat(top1[i],geneSymbol1[i],GOAnot,substr(Term(GOAnot),1,30), "\n")
  }
}

## 1450826_a_at Saa3 GO:0006953 acute-phase response 
## 1450826_a_at Saa3 GO:0007252 I-kappaB phosphorylation 
## 1450826_a_at Saa3 GO:0009617 response to bacterium 
## 1450826_a_at Saa3 GO:0035634 response to stilbenoid 
## 1450826_a_at Saa3 GO:0060326 cell chemotaxis 
## 1450826_a_at Saa3 GO:0071347 cellular response to interleuk 
## 1450826_a_at Saa3 GO:0005576 extracellular region 
## 1450826_a_at Saa3 GO:0005615 extracellular space 
## 1450826_a_at Saa3 GO:0034364 high-density lipoprotein parti 
## 1450826_a_at Saa3 GO:0035662 Toll-like receptor 4 binding 
## 1450826_a_at Saa3 GO:0042056 chemoattractant activity

Como se ve en el ejemplo el número de anotaciones para un gen en la Gene Ontology es muy alto, aparte de que, aunque aquí no se muestra, no todas las anotaciones tienen la misma fiabilidad.

Aparte del problema de lo extenso de las anotaciones está el hecho de que antes de empezar a hacer inferencias sobre el significado de una anotación debería poderse establecer si dicha anotación está relacionada con el proceso que se está estudiando o aparece por azar entre la muchas anotaciones de los genes de la lista. Hay diferentes métodos y modelos para hacer esto (ver Draghici y colegas, Khatri and Drăghici (2005) o Mosquera and Sánchez–Pla, Mosquera and Sánchez-Pla (2005),Sanchez:2007b) pero aquí se presentará brevemente lo que se conoce por Análisis de enriquecimiento.

8.5.1 Análisis de enriquecimiento

El objetivo del análisis de enriquecimiento es establecer si una determinada categoría, que representa, por ejemplo, un proceso biológico (GO) o una vía (KEGG), aparece más (“enriquecido”) o menos (“pobre”) a menudo en la lista de genes seleccionados que en la población (génica), desde donde se han obtenido, es decir, el array, el genoma, o simplemente los genes que fueron seleccionados para la prueba. La idea básica es que si una función aparece más a menudo en la lista que en general es probable que tenga algo que ver con el proceso en base al cual se ha seleccionado la lista que se estudia.

Por ejemplo, consideremos un experimento que da una lista de genes y el 10% de los genes más diferencialmente expresados están asociados con el término apoptosis en la GO (GO:0006915). Esto puede parecer una proporción inusualmente grande de la lista de genes, dado que la apoptosis es un proceso biológico muy específico. Para determinar cuánto más grande de lo normal es esta proporción, debe ser comparada con la proporción de genes relacionados con apoptosis en la lista de genes de referencia, que suele ser el conjunto de todos los genes del microarray.

El análisis estadístico realizado para comparar proporciones es un test Hipergeométrico o el test exacto de Fisher que se utiliza para probar la hipótesis:

\[ H_0: \ p_{sel}^A = p_{all}^A \ vs \ H_1: \ p_{sel}^A \neq p_{all}^A, \]

donde $A$ representa el conjunto de genes cuya más/menos representación está siendo considerada, $p_{sel}^A$ es la proporción de genes seleccionados que están incluidos en este conjunto de genes y $p_{all}^A$ es la proporción de genes de la lista de referencia.

8.5.2 Análisis de enriquecimiento con Bioconductor

Hay muchas herramientas que pueden ayudar a realizar este análisis. Siguiendo nuestra costumbre usaremos las que contiene el paquete de R(Bioconductor), principalmente en el paquete GOstats.

El análisis realizado utilizando este paquete procede de la forma siguiente:

Toma como entrada los identificadores de Entrez o de Affy de la lista de genes seleccionada así como el nombre del paquete de la anotación correspondiente al array que ha sido usado para el análisis.
La salida del análisis es la lista de categorias que aparece más/menos representado en cada conjunto seleccionado.

El código siguiente ilustra como se procederá para realizar un análisis de enriquecimiento con las listas de genes seleccionados en la primera de las comparaciones realizadas en este capítulo.

require(GOstats)
require(mouse4302.db)
require(org.Mm.eg.db)

  # Seleccionamos la "topTable"
  topTab <- topTab_LPS.in.AGED
  # Definimos el universo de genes: todos los que se han incluido en el análisis
  # EL programa trabaja con identificadores "entrez" y no admite duplicados

  entrezUniverse <- unique(getEG(as.character(rownames(topTab)), "mouse4302.db"))

  # Escogemos los grupos de sondas a incluir en el análisis
  # Este análisis trabaja bien con varios centenares de genes
  # por lo que es habitual basarse en p-valores sin ajustar para incluirlos

  whichGenes<-topTab["adj.P.Val"]<0.001
  geneIds <-   unique(getEG(as.character(rownames(topTab)[whichGenes]),"mouse4302.db"))

  # Creamos los "hiperparámetros" en que se basa el análisis
  GOparams = new("GOHyperGParams",
    geneIds=geneIds, universeGeneIds=entrezUniverse,
    annotation="org.Mm.eg.db", ontology="BP",
    pvalueCutoff=0.001, conditional=FALSE,
    testDirection="over")

  # Ejecutamos los análisis

  GOhyper = hyperGTest(GOparams)

# Creamos un informe html con los resultados
   comparison = "topTab_LPS.in.AGED"
   GOfilename =file.path(resultsDir,
     paste("GOResults.",comparison,".html", sep=""))
  htmlReport(GOhyper, file = GOfilename, summary.args=list("htmlLinks"=TRUE))

El análisis basado en la base de datos de Pathways, KEGG será básicamente el mismo, cambiando únicamente el tipo de “hiperparámetro” invocado.

  KEGGparams = new("KEGGHyperGParams",
    geneIds=geneIds, universeGeneIds=entrezUniverse,
    annotation="org.Mm.eg.db",
    pvalueCutoff=0.01, testDirection="over")

  # Ejecutamos los análisis

  KEGGhyper = hyperGTest(KEGGparams)

# Creamos un informe html con los resultados
 comparison = "topTab_LPS.in.AGED"
 KEGGfilename =file.path(resultsDir,
     paste("KEGGResults.",comparison,".html", sep=""))
  htmlReport(KEGGhyper, file = KEGGfilename, summary.args=list("htmlLinks"=TRUE))

9 Referencias

Alizadeh, A., M. B. Eisen, E. Davis, C. Ma, I. Lossos, A. Rosenwald, J. Boldrick, et al. 2000. “Distinct Types of Diffuse Large B–Cell Lymphoma Identified by Gene Expression Profiling.” Nature 403: 503–11.

Allison, David B, Xiangqin Cui, Grier P Page, and Mahyar Sabripour. 2006. “Microarray Data Analysis: From Disarray to Consolidation and Consensus.” Nat Rev Genet 7 (1): 55–65.

Barrier, Alain, Pierre-Yves Boelle, Antoinette Lemoine, Antoine Flahault, Sandrine Dudoit, and Michel Huguier. 2007. “[Gene Expression Profiling in Colon Cancer].” Bull Acad Natl Med 191 (6): 1091-101; discussion 1102-3.

Benjamini, Y., and Y. Hochberg. 1995. “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” Journal of the Royal Statistical Society B 57: 289–300.

Callow, M. J., S. Dudoit, E. L. Gong, T. P. Speed, and E. M. Rubin. 2000. “Microarray Expression Profiling Identifies Genes with Altered Expression in HDL-Deficient Mice.” Genome Research.

Chelvarajan, R. L., Y. Liu, D. Popa, M. L. Getchell, T. V. Getchell, A. J. Stromberg, and S. Bondada. 2006. “Molecular basis of age-associated cytokine dysregulation in LPS-stimulated macrophages.” Journal of Leukocyte Biology 79 (6): 1314.

Dudoit, S., J. P. Shaffer, and J. C. Boldrick. 2003. “Multiple Hypothesis Testing in Microarray Experiments.” Statistical Science 18: 71–103.

Dyrskjøt, L., T. Thykjaer, M. Kruhøffer, J. L. Jensen, N. Marcussen, S. Hamilton-Dutoit, H. Wolf, and T. F. Ørntoft. 2002. “Identifying distinct classes of bladder carcinoma using microarrays.” Nature Genetics 33: 90–96.

Faraway, Julian J. 2004. Linear Models with r. 1st ed. Chapman; Hall/CRC.

Geschwind, D. H., and J. P. Gregg. 2002. Microarrays for the neurosciences: an essential guide. MIT Press.

Golub, T. R., D. K. Slonim, P. Tamayo, C. Huard, M. Gaasenbeek, J. P. Mesirov, H. Coller, et al. 1999. “Molecular Classification of Cancer: Class Discovery and Class Prediction by Gene Expression Monitoring.” Science 286: 531–37.

Hedenfalk, I. A., M. Ringnér, J. M. Trent, and A. Borg. 2002. “Gene Expression in Inherited Breast Cancer.” Adv Cancer Res 84: 1–34. http://view.ncbi.nlm.nih.gov/pubmed/11883525.

Imbeaud, Sandrine, and Charles Auffray. 2005. “’The 39 Steps’ in Gene Expression Profiling: Critical Issues and Proposed Best Practices for Microarray Experiments.” Drug Discovery Today 10 (17): 1175–82. https://doi.org/10.1016/S1359-6446(05)03565-8.

Kerr, M K, and G A Churchill. 2001. “Experimental Design for Gene Expression Microarrays.” Biostatistics.

Khatri, P., and S. Drăghici. 2005. “Ontological Analysis of Gene Expression Data: Current Tools, Limitations, and Problems.” Bioinformatics 18: 3587–95.

Kupin, Isabelle Lesur. 2005. “Study of the Transcriptome of the Prematurely Aging Dna-2 Yeast Mutant Using a New System Allowing Comparative DNA Microarray Analysis.” PhD thesis, Universite Bordeaux I.

Lee, M. L. T., and GA Whitmore. 2002. “Power and sample size for DNA microarray studies.” Statistics in Medicine 21 (23): 3543–70.

Lin, Simon M, Jyothi Devakumar, and Warren A Kibbe. 2006. “Improved Prediction of Treatment Response Using Microarrays and Existing Biological Knowledge.” Pharmacogenomics 7 (3): 495–501. https://doi.org/10.2217/14622416.7.3.495.

Moore, Shanna, Julia Vrebalov, Paxton Payton, and Jim Giovannoni. 2002. “Use of Genomics Tools to Isolate Key Ripening Genes and Analyse Fruit Maturation in Tomato.” Journal of Experimental Botany 53 (377): 2023–30. https://doi.org/10.1093/jxb/erf057.

Mosquera, J-L., and A. Sánchez-Pla. 2005. “A Comparative Study of GO Mining Programs.” In X Conferencia Española de Biometría. Sociedad Española de Biometría.

Quackenbush, J. 2002. “Microarray Data Normalization and Transformation.” Nature Genet. 32: 496–501.

Schena, M. 1999. Microarray Analysis. J. Wiley & Sons, Hoboken, NJ.

Scholtens, Denise, Alexander Miron, Faisal M. Merchant, Arden Miller, Penelope L. Miron, J. Dirk Iglehart, and Robert Gentleman. 2004. “Analyzing Factorial Designed Microarray Experiments.” J. Multivar. Anal. 90 (1): 19–43. https://doi.org/http://dx.doi.org/10.1016/j.jmva.2004.02.004.

Simon, Richard M., Edward L. Korn, Lisa M. McShane, Michael D. Radmacher, George W. Wright, and Yingdong Zhao. 2003. Design and Analysis of DNA Microarray Investigations. Springer-Verlag.

Smyth, Gordon K, Joëlle Michaud, and Hamish S Scott. 2005. “Use of Within-Array Replicate Spots for Assessing Differential Expression in Microarray Experiments.” Bioinformatics 21 (9): 2067–75.

Speed, T. 2003. Statistical Analysis of Gene Expression Data. Boca Raton, Fla.: Chapman & Hall/CRC.

Tibshirani, R. 2006. “A simple method for assessing sample sizes in microarray experiments.” BMC Bioinformatics 7 (1): 106.

van ’t Veer, Laura J, Hongyue Dai, Marc J van de Vijver, Yudong D He, Augustinus A M Hart, Mao Mao, Hans L Peterse, et al. 2002. “Gene Expression Profiling Predicts Clinical Outcome of Breast Cancer.” Nature 415 (6871): 530–36.

Wang, Zhong, Mark Gerstein, and Michael Snyder. 2009. “RNA-Seq: A Revolutionary Tool for Transcriptomics.” Nat Rev Genet 10 (1): 57–63. https://doi.org/10.1038/nrg2484.

Yang, Y. H., M. J. Buckley, S. Dudoit, and T. P. Speed. 2002. “Comparison of Methods for Image Analysis on cDNA Microarray Data.” Journal of Computational and Graphical Statistic 11 (1).

Zhu, Qianqian, Jeffrey Miecznikowski, and Marc Halfon. 2010. “Preferred Analysis Methods for Affymetrix GeneChips. II. An Expanded, Balanced, Wholly-Defined Spike-in Dataset.” BMC Bioinformatics 11 (1): 285. https://doi.org/10.1186/1471-2105-11-285.

Analisis de datos ómicos

Alex Sánchez-Pla\(^{1,2}\) y Ricardo Gonzalo Sanz\(^2\)

\(^1\) Departament de Genètica, Microbiologia i Estadística. Universitat de Barcelona (UB).

\(^2\) Unitat de Estadística i Bioinformàtica. Vall d’Hebron Institut de Recerca (VHIR).

1 Presentación

1.1 Instalación de paquetes de R y código de los capítulos

2 Introduccion a los microarrays y otras tecnologías omicas

2.1 Antecedentes históricos

2.2 Los microarrays

2.2.1 Aplicaciones de los microarrays

2.3 Cómo funcionan los microarrays

2.3.1 Microarrays de dos colores

2.3.2 Microarrays de un color

2.3.2.1 Fabricación de los arrays de oligonucleótidos

2.3.3 Realización de un experimento de microarrays

2.3.4 Como se mide la expresión

2.3.4.1 Medición de la expresión relativa en arrays de dos colores

2.3.4.2 Medición de la expresión absoluta en arrays de un color

2.4 Bioinformática de Microarrays

2.4.1 Software para el análisis de datos de microarrays

2.4.1.1 Software libre o de código abierto

2.4.1.2 El proyecto bioconductor

2.4.1.3 El software propietario

2.4.2 Bases de datos de Microarrays

2.5 Extensiones (1): Otros tipos de microarrays

2.5.1 Otros microarrays de ADN

2.5.1.1 Arrays de genotipado (SNPs)

2.5.2 Otros tipos de microarrays: Proteínas o carbohidratos

2.6 Extensiones (2): NGS: Next(Now) Generation Sequencing

2.6.1 Visión global de las tecnologías de secuenciación

2.6.2 De la secuenciación Sanger a la ultrasecuenciación

2.6.3 Tecnologías de ultrasecuenciación

3 El proceso de análisis de microarrays

3.1 Tipos de estudios

3.1.1 Comparación de grupos o “Class comparison”

3.1.2 Predicción de clase o “Class prediction”

3.1.3 Descubrimiento de clases o “Class discovery”

3.1.4 Otros tipos de estudios

3.2 Algunos ejemplos concretos

3.2.1 Estudio de procesos regulados por citoquinas

3.2.2 Clasificación molecular de la leucemia

3.2.3 Efecto del estrogeno y el tiempo de administración

3.2.4 Efecto del CCL4 en la expresión génica

3.2.5 Estudio de la mutación “swirl” en el desarrollo del pez cebra

3.2.6 Análisis de patrones en el ciclo celular

3.2.7 Recapitulación

3.3 El proceso de análisis de microarrays

4 Diseño de experimentos de microarrays

4.1 Fuentes de variabilidad

4.2 Tipos de variabilidad

4.3 Conceptos básicos de diseño de Experimentos

4.4 Principios de diseño de experimentos

4.4.1 Replicación

4.4.1.1 Potencia y tamaño de la muestra

4.4.1.2 Pooling

4.4.2 Aleatorización

4.4.3 Bloqueo o control local

4.5 Diseños experimentales para microarrays de dos colores

5 Exploración de los datos, control de calidad y preprocesado

5.1 Introducción

5.1.1 Nivel de análisis y tipo de microarray

5.1.2 Datos de partida

5.1.2.1 Arrays de dos colores (o de “cDNA”)

5.1.2.2 Arrays de un color (Affymetrix)

5.1.2.3 Datos de ejemplo

5.2 Exploración y control de calidad

5.2.1 Control de calidad con gráficos estadísticos generales

5.2.1.1 Histogramas y gráficos de densidad

5.2.1.2 Diagramas de caja o “boxplots”

5.2.1.3 Gráficos de componentes principales

5.2.1.4 Imagen del chip

5.2.2 Gráficos de diagnóstico para microarrays de dos colores

5.2.2.1 Diagramas de dispersión y “MA-plots”

5.2.2.2 Imagen del array

5.2.2.3 Histogramas de señales y de la relación señal–ruído

5.2.2.4 Boxplots

5.2.3 Gráficos de diagnóstico para microarrays de un color

5.2.3.1 Imagen del chip

5.2.3.2 Gráfico “M-A”

5.2.3.3 Modelos de bajo nivel (“Probe-Level-Models” o PLM)

6.5 Modelos lineales para la selección de genes: `limma`

8.2.1 Ejemplo: comparación de listas en el caso `celltypes`